From 48f2b5d59b0dac8882e1bba9946f9d9c308eccbe Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Sun, 5 May 2024 16:28:08 +0800 Subject: [PATCH] Refactor: optimized cubic spline (#4081) * redesign eval to enable potential vectorization * fix interface * refactored version init commit * clean up * clean up * fix & a few tests * theoretical error bound test * clean up * better unit test for cross-check * systematic cross-check * minor change * clean up * clean up * minor internal change * add explanation --------- Co-authored-by: Mohan Chen --- source/module_base/cubic_spline.cpp | 711 ++++--- source/module_base/cubic_spline.h | 749 ++++--- source/module_base/test/cubic_spline_test.cpp | 1755 ++++------------- source/module_base/test/data/cos_periodic.dat | 7 + .../module_base/test/data/exp_first_deriv.dat | 7 + source/module_base/test/data/gen_ref.py | 64 + .../test/data/log_second_deriv.dat | 7 + .../module_base/test/data/sin_not_a_knot.dat | 7 + source/module_base/test/data/sqrt_mix_bc.dat | 7 + .../test/data/three_points_not_a_knot.dat | 7 + .../test/data/two_points_first_deriv.dat | 7 + .../test/data/two_points_periodic.dat | 7 + .../test/data/two_points_second_deriv.dat | 7 + .../module_nao/numerical_radial.cpp | 3 +- source/module_basis/module_nao/projgen.cpp | 6 +- .../module_nao/two_center_table.cpp | 8 +- 16 files changed, 1487 insertions(+), 1872 deletions(-) create mode 100644 source/module_base/test/data/cos_periodic.dat create mode 100644 source/module_base/test/data/exp_first_deriv.dat create mode 100644 source/module_base/test/data/gen_ref.py create mode 100644 source/module_base/test/data/log_second_deriv.dat create mode 100644 source/module_base/test/data/sin_not_a_knot.dat create mode 100644 source/module_base/test/data/sqrt_mix_bc.dat create mode 100644 source/module_base/test/data/three_points_not_a_knot.dat create mode 100644 source/module_base/test/data/two_points_first_deriv.dat create mode 100644 source/module_base/test/data/two_points_periodic.dat create mode 100644 source/module_base/test/data/two_points_second_deriv.dat diff --git a/source/module_base/cubic_spline.cpp b/source/module_base/cubic_spline.cpp index 1bab615658..efe58ee694 100644 --- a/source/module_base/cubic_spline.cpp +++ b/source/module_base/cubic_spline.cpp @@ -1,357 +1,554 @@ #include "cubic_spline.h" -#include #include -#include -#include +#include +#include #include #include -namespace ModuleBase +using ModuleBase::CubicSpline; + +extern "C" { + // solve a tridiagonal linear system + void dgtsv_(int* N, int* NRHS, double* DL, double* D, double* DU, double* B, int* LDB, int* INFO); +}; -CubicSpline::CubicSpline(CubicSpline const& other) + +CubicSpline::BoundaryCondition::BoundaryCondition(BoundaryType type) + : type(type) { - n_ = other.n_; - is_uniform_ = other.is_uniform_; - uniform_thr_ = other.uniform_thr_; + assert(type == BoundaryType::periodic || type == BoundaryType::not_a_knot); +} - x_ = new double[n_]; - y_ = new double[n_]; - s_ = new double[n_]; - std::memcpy(x_, other.x_, n_ * sizeof(double)); - std::memcpy(y_, other.y_, n_ * sizeof(double)); - std::memcpy(s_, other.s_, n_ * sizeof(double)); +CubicSpline::BoundaryCondition::BoundaryCondition(BoundaryType type, double val) + : type(type), val(val) +{ + assert(type == BoundaryType::first_deriv || type == BoundaryType::second_deriv); } -CubicSpline& CubicSpline::operator=(CubicSpline const& other) + +CubicSpline::CubicSpline( + int n, + const double* x, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end +): n_spline_(1), n_(n), xmin_(x[0]), xmax_(x[n - 1]), x_(x, x + n), y_(2 * n) { - if (this != &other) + std::copy(y, y + n, y_.begin()); + build(n, x, y, bc_start, bc_end, &y_[n]); +} + + +CubicSpline::CubicSpline( + int n, + double x0, + double dx, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end +): n_spline_(1), n_(n), xmin_(x0), xmax_(x0 + (n - 1) * dx), dx_(dx), y_(2 * n) +{ + std::copy(y, y + n, y_.begin()); + build(n, dx, y, bc_start, bc_end, &y_[n]); +} + + +void CubicSpline::add( + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end +) +{ + int offset = n_spline_ * 2 * n_; + y_.resize(offset + 2 * n_); + + std::copy(y, y + n_, &y_[offset]); + double* dy = &y_[offset + n_]; + if (x_.empty()) // evenly spaced knots { - cleanup(); - n_ = other.n_; - is_uniform_ = other.is_uniform_; - uniform_thr_ = other.uniform_thr_; - - x_ = new double[n_]; - y_ = new double[n_]; - s_ = new double[n_]; - - std::memcpy(x_, other.x_, n_ * sizeof(double)); - std::memcpy(y_, other.y_, n_ * sizeof(double)); - std::memcpy(s_, other.s_, n_ * sizeof(double)); + build(n_, dx_, y, bc_start, bc_end, dy); } + else + { + build(n_, x_.data(), y, bc_start, bc_end, dy); + } + ++n_spline_; +} + - return *this; +void CubicSpline::eval( + int n_interp, + const double* x_interp, + double* y_interp, + double* dy_interp, + double* d2y_interp, + int i_spline +) const +{ + assert(0 <= i_spline && i_spline < n_spline_); + + const double* y = &y_[i_spline * 2 * n_]; + const double* dy = y + n_; + if (x_.empty()) // evenly spaced knots + { + eval(n_, xmin_, dx_, y, dy, n_interp, x_interp, y_interp, dy_interp, d2y_interp); + } + else + { + eval(n_, x_.data(), y, dy, n_interp, x_interp, y_interp, dy_interp, d2y_interp); + } } -void CubicSpline::cleanup() + +void CubicSpline::multi_eval( + int n_spline, + const int* i_spline, + double x_interp, + double* y_interp, + double* dy_interp, + double* d2y_interp +) const { - delete[] x_; - delete[] y_; - delete[] s_; + assert(std::all_of(i_spline, i_spline + n_spline, + [this](int i) { return 0 <= i && i < n_spline_; })); + _validate_eval(n_, {xmin_, dx_}, x_.empty() ? nullptr : x_.data(), + y_.data(), &y_[n_], 1, &x_interp); + + int p = 0; + double dx = 0.0, r = 0.0; + if (x_.empty()) // evenly spaced knots + { + p = _index(n_, xmin_, dx_, x_interp); + dx = dx_; + r = (x_interp - xmin_) / dx - p; + } + else + { + p = _index(n_, x_.data(), x_interp); + dx = x_[p + 1] - x_[p]; + r = (x_interp - x_[p]) / dx; + } + + const double r2 = r * r; + const double r3 = r2 * r; + double wy0 = 0.0, wy1 = 0.0, ws0 = 0.0, ws1 = 0.0; + int offset = 0; + + if (y_interp) + { + wy1 = 3.0 * r2 - 2.0 * r3; + wy0 = 1.0 - wy1; + ws0 = (r - 2.0 * r2 + r3) * dx; + ws1 = (r3 - r2) * dx; + for (int i = 0; i < n_spline; ++i) + { + offset = i_spline[i] * 2 * n_ + p; + y_interp[i] = wy0 * y_[offset] + wy1 * y_[offset + 1] + + ws0 * y_[offset + n_] + ws1 * y_[offset + n_ + 1]; + } + } + + if (dy_interp) + { + wy1 = 6.0 * (r - r2) / dx; // wy0 = -wy1 + ws0 = 3.0 * r2 - 4.0 * r + 1.0; + ws1 = 3.0 * r2 - 2.0 * r; + for (int i = 0; i < n_spline; ++i) + { + offset = i_spline[i] * 2 * n_ + p; + dy_interp[i] = wy1 * (y_[offset + 1] - y_[offset]) + + ws0 * y_[offset + n_] + ws1 * y_[offset + n_ + 1]; + } + } + + if (d2y_interp) + { + wy1 = (6.0 - 12.0 * r) / (dx * dx); // wy0 = -wy1 + ws0 = (6.0 * r - 4.0) / dx; + ws1 = (6.0 * r - 2.0) / dx; + for (int i = 0; i < n_spline; ++i) + { + offset = i_spline[i] * 2 * n_ + p; + d2y_interp[i] = wy1 * (y_[offset + 1] - y_[offset]) + + ws0 * y_[offset + n_] + ws1 * y_[offset + n_ + 1]; + } + } +} - x_ = nullptr; - y_ = nullptr; - s_ = nullptr; + +void CubicSpline::multi_eval( + double x, + double* y, + double* dy, + double* d2y +) const +{ + std::vector i_spline(n_spline_); + std::iota(i_spline.begin(), i_spline.end(), 0); + multi_eval(i_spline.size(), i_spline.data(), x, y, dy, d2y); } -void CubicSpline::check_build(const int n, - const double* const x, - const double* const y, - BoundaryCondition bc_start, - BoundaryCondition bc_end) + +void CubicSpline::build( + int n, + const double* x, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end, + double* dy +) { - assert(n > 1); + std::vector dx(n); + std::adjacent_difference(x, x + n, dx.begin()); + _build(n, &dx[1], y, bc_start, bc_end, dy); +} - // periodic boundary condition must apply to both ends - assert((bc_start == BoundaryCondition::periodic) == (bc_end == BoundaryCondition::periodic)); - // y[0] must equal y[n-1] for periodic boundary condition - assert(bc_start != BoundaryCondition::periodic || y[0] == y[n - 1]); +void CubicSpline::build( + int n, + double dx, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end, + double* dy +) +{ + std::vector dx_(n - 1, dx); + _build(n, dx_.data(), y, bc_start, bc_end, dy); +} + + +void CubicSpline::eval( + int n, + const double* x, + const double* y, + const double* dy, + int n_interp, + const double* x_interp, + double* y_interp, + double* dy_interp, + double* d2y_interp +) +{ + _validate_eval(n, {}, x, y, dy, n_interp, x_interp); + + // indices of the polynomial segments that contain x_interp + std::vector _ind(n_interp); + std::transform(x_interp, x_interp + n_interp, _ind.begin(), + [n, x](double x_i) { return _index(n, x, x_i); }); - // if one of the boundary condition is not-a-knot, n must be larger than 2 - assert((bc_start != BoundaryCondition::not_a_knot && bc_end != BoundaryCondition::not_a_knot) || n > 2); + std::vector buffer(n_interp * 5); + double* _w = buffer.data(); + double* _c0 = _w + n_interp; + double* _c1 = _c0 + n_interp; + double* _c2 = _c1 + n_interp; + double* _c3 = _c2 + n_interp; - for (int i = 0; i != n - 1; ++i) + for (int i = 0; i < n_interp; ++i) { - assert(x[i + 1] > x[i] && "CubicSpline: x must be strictly increasing"); + int p = _ind[i]; + double dx = x[p + 1] - x[p]; + double inv_dx = 1.0 / dx; + double dd = (y[p + 1] - y[p]) * inv_dx; + _w[i] = x_interp[i] - x[p]; + _c0[i] = y[p]; + _c1[i] = dy[p]; + _c3[i] = (_c1[i] + dy[p + 1] - 2.0 * dd) * inv_dx * inv_dx; + _c2[i] = (dd - _c1[i]) * inv_dx - _c3[i] * dx; } + + _cubic(n_interp, _w, _c0, _c1, _c2, _c3, y_interp, dy_interp, d2y_interp); } -void CubicSpline::check_interp(const int n, - const double* const x, - const double* const y, - const double* const s, - const int n_interp, - const double* const x_interp, - double* const y_interp, - double* const dy_interp) + +void CubicSpline::eval( + int n, + double x0, + double dx, + const double* y, + const double* dy, + int n_interp, + const double* x_interp, + double* y_interp, + double* dy_interp, + double* d2y_interp +) +{ + _validate_eval(n, {x0, dx}, nullptr, y, dy, n_interp, x_interp); + + // indices of the polynomial segments that contain x_interp + std::vector _ind(n_interp); + std::transform(x_interp, x_interp + n_interp, _ind.begin(), + [n, x0, dx](double x_i) { return _index(n, x0, dx, x_i); }); + + std::vector buffer(n_interp * 5); + double* _w = buffer.data(); + double* _c0 = _w + n_interp; + double* _c1 = _c0 + n_interp; + double* _c2 = _c1 + n_interp; + double* _c3 = _c2 + n_interp; + + double inv_dx = 1.0 / dx; + double inv_dx2 = inv_dx * inv_dx; + for (int i = 0; i < n_interp; ++i) + { + int p = _ind[i]; + double dd = (y[p + 1] - y[p]) * inv_dx; + _w[i] = x_interp[i] - x0 - p * dx; + _c0[i] = y[p]; + _c1[i] = dy[p]; + _c3[i] = (_c1[i] + dy[p + 1] - 2.0 * dd) * inv_dx2; + _c2[i] = (dd - _c1[i]) * inv_dx - _c3[i] * dx; + } + + _cubic(n_interp, _w, _c0, _c1, _c2, _c3, y_interp, dy_interp, d2y_interp); +} + + +void CubicSpline::_validate_build( + int n, + const double* dx, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end +) +{ + assert(n > 1); + + // if periodic boundary condition is specified, it must be applied to both ends + assert((bc_start.type == BoundaryType::periodic) + == (bc_end.type == BoundaryType::periodic)); + + // y[0] must equal y[n-1] for periodic boundary condition + assert(bc_start.type != BoundaryType::periodic || y[0] == y[n - 1]); + + // not-a-knot boundary condition requires the existence of "internal" knot + // so n must be at least 3 + assert((bc_start.type != BoundaryType::not_a_knot && + bc_end.type != BoundaryType::not_a_knot) || n > 2); + + // knots must be strictly increasing + assert(std::all_of(dx, dx + n - 1, [](double d) { return d > 0.0; })); +} + + +void CubicSpline::_validate_eval( + int n, + const double (&u)[2], + const double* x, + const double* y, + const double* dy, + int n_interp, + const double* x_interp +) { - assert(n > 1 && x && y && s); // make sure the interpolant exists - assert(n_interp > 0 && x_interp); // make sure the interpolation points exist - assert(y_interp || dy_interp); // make sure at least one of y or dy is not null + assert(n > 1 && y && dy); + assert((x && std::is_sorted(x, x + n, std::less_equal())) || u[1] > 0.0); - // check that x_interp is in the range of the interpolant + assert((n_interp > 0 && x_interp) || n_interp == 0); + + double xmin = x ? x[0] : u[0]; + double xmax = x ? x[n - 1] : u[0] + (n - 1) * u[1]; assert(std::all_of(x_interp, x_interp + n_interp, - [n, x](const double xi) { return xi >= x[0] && xi <= x[n - 1]; })); + [xmin, xmax](double x_i) { return xmin <= x_i && x_i <= xmax; })); } -void CubicSpline::build(const int n, - const double* const x, - const double* const y, - double* const s, - BoundaryCondition bc_start, - BoundaryCondition bc_end, - const double deriv_start, - const double deriv_end) + +void CubicSpline::_build( + int n, + const double* dx, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end, + double* dy +) { - check_build(n, x, y, bc_start, bc_end); + _validate_build(n, dx, y, bc_start, bc_end); - if (n == 2 && bc_start == BoundaryCondition::periodic) - { // in this case the polynomial is a constant - s[0] = s[1] = 0.0; + if (n == 2 && bc_start.type == BoundaryType::periodic) + { + dy[0] = dy[1] = 0.0; // the only possible solution: constant } - else if (n == 3 && bc_start == BoundaryCondition::not_a_knot && bc_end == BoundaryCondition::not_a_knot) - { // in this case two conditions coincide; simply build a parabola that passes through the three data points - double idx10 = 1. / (x[1] - x[0]); - double idx21 = 1. / (x[2] - x[1]); - double idx20 = 1. / (x[2] - x[0]); - - s[0] = -y[0] * (idx10 + idx20) + y[1] * (idx21 + idx10) + y[2] * (idx20 - idx21); - s[1] = -y[1] * (-idx10 + idx21) + y[0] * (idx20 - idx10) + y[2] * (idx21 - idx20); - s[2] = s[1] + 2.0 * (-y[1] * idx10 + y[2] * idx20) + 2.0 * y[0] * idx10 * idx20 * (x[2] - x[1]); + else if (n == 3 && bc_start.type == BoundaryType::not_a_knot + && bc_end.type == BoundaryType::not_a_knot) + { + // in this case two conditions coincide + // simply build a parabola that passes through the three data points + double dd01 = (y[1] - y[0]) / dx[0]; // divided difference f[x0,x1] + double dd12 = (y[2] - y[1]) / dx[1]; // f[x1,x2] + double dd012 = (dd12 - dd01) / (dx[0] + dx[1]); // f[x0,x1,x2] + + dy[0] = dd01 - dd012 * dx[0]; + dy[1] = 2.0 * dd01 - dy[0]; + dy[2] = dd01 + dd012 * (dx[0] + 2.0 * dx[1]); } else { - // double* dx = new double[n - 1]; - // double* dd = new double[n - 1]; - std::vector dx(n - 1); - std::vector dd(n - 1); - - // tridiagonal linear system (cyclic if using the periodic boundary condition) - // double* diag = new double[n]; - // double* subdiag = new double[n - 1]; - // double* supdiag = new double[n - 1]; - std::vector diag(n); - std::vector subdiag(n - 1); - std::vector supdiag(n - 1); - - for (int i = 0; i != n - 1; ++i) - { - dx[i] = x[i + 1] - x[i]; - dd[i] = (y[i + 1] - y[i]) / dx[i]; - } + std::vector buffer(4 * n); + double* dd = buffer.data(); // divided differences + std::adjacent_difference(y, y + n, dd); + dd += 1; // the first element computed by adjacent_difference is not a difference + std::transform(dd, dd + n - 1, dx, dd, std::divides()); + + // tridiagonal linear system (cyclic tridiagonal if periodic boundary condition) + double* d = buffer.data() + n; // main diagonal + double* l = buffer.data() + 2 * n; // subdiagonal + double* u = buffer.data() + 3 * n; // superdiagonal + + //*********************************************** // common part of the tridiagonal linear system + //*********************************************** + std::copy(dx + 1, dx + n - 1, l); + std::copy(dx, dx + n - 2, u + 1); + for (int i = 1; i != n - 1; ++i) { - diag[i] = 2.0 * (dx[i - 1] + dx[i]); - supdiag[i] = dx[i - 1]; - subdiag[i - 1] = dx[i]; - s[i] = 3.0 * (dd[i - 1] * dx[i] + dd[i] * dx[i - 1]); + d[i] = 2.0 * (dx[i - 1] + dx[i]); + dy[i] = 3.0 * (dd[i - 1] * dx[i] + dd[i] * dx[i - 1]); } - if (bc_start == BoundaryCondition::periodic) + //*********************************************** + // boundary-specific part + //*********************************************** + if (bc_start.type == BoundaryType::periodic) { // exclude s[n-1] and solve a a cyclic tridiagonal linear system of size n-1 - diag[0] = 2.0 * (dx[n - 2] + dx[0]); - supdiag[0] = dx[n - 2]; - subdiag[n - 2] = dx[0]; - s[0] = 3.0 * (dd[0] * dx[n - 2] + dd[n - 2] * dx[0]); - solve_cyctri(n - 1, diag.data(), supdiag.data(), subdiag.data(), s); - s[n - 1] = s[0]; + d[0] = 2.0 * (dx[n - 2] + dx[0]); + u[0] = dx[n - 2]; + l[n - 2] = dx[0]; + dy[0] = 3.0 * (dd[0] * dx[n - 2] + dd[n - 2] * dx[0]); + _solve_cyctri(n - 1, d, u, l, dy); + dy[n - 1] = dy[0]; } else { - switch (bc_start) + switch (bc_start.type) { - case BoundaryCondition::first_deriv: - diag[0] = 1.0 * dx[0]; - supdiag[0] = 0.0; - s[0] = deriv_start * dx[0]; + case BoundaryType::first_deriv: + d[0] = 1.0 * dx[0]; + u[0] = 0.0; + dy[0] = bc_start.val * dx[0]; break; - case BoundaryCondition::second_deriv: - diag[0] = 2.0 * dx[0]; - supdiag[0] = 1.0 * dx[0]; - s[0] = (3.0 * dd[0] - 0.5 * deriv_start * dx[0]) * dx[0]; + case BoundaryType::second_deriv: + d[0] = 2.0 * dx[0]; + u[0] = 1.0 * dx[0]; + dy[0] = (3.0 * dd[0] - 0.5 * bc_start.val * dx[0]) * dx[0]; break; default: // BoundaryCondition::not_a_knot - diag[0] = dx[1]; - supdiag[0] = x[2] - x[0]; - s[0] = (dd[0] * dx[1] * (dx[0] + 2 * (x[2] - x[0])) + dd[1] * dx[0] * dx[0]) / (x[2] - x[0]); + d[0] = dx[1]; + u[0] = dx[0] + dx[1]; + dy[0] = (dd[0] * dx[1] * (dx[0] + 2 * u[0]) + dd[1] * dx[0] * dx[0]) / u[0]; } - switch (bc_end) + switch (bc_end.type) { - case BoundaryCondition::first_deriv: - diag[n - 1] = 1.0 * dx[n - 2]; - subdiag[n - 2] = 0.0; - s[n - 1] = deriv_end * dx[n - 2]; + case BoundaryType::first_deriv: + d[n - 1] = 1.0 * dx[n - 2]; + l[n - 2] = 0.0; + dy[n - 1] = bc_end.val * dx[n - 2]; break; - case BoundaryCondition::second_deriv: - diag[n - 1] = 2.0 * dx[n - 2]; - subdiag[n - 2] = 1.0 * dx[n - 2]; - s[n - 1] = (3.0 * dd[n - 2] + 0.5 * deriv_end * dx[n - 2]) * dx[n - 2]; + case BoundaryType::second_deriv: + d[n - 1] = 2.0 * dx[n - 2]; + l[n - 2] = 1.0 * dx[n - 2]; + dy[n - 1] = (3.0 * dd[n - 2] + 0.5 * bc_end.val * dx[n - 2]) * dx[n - 2]; break; default: // BoundaryCondition::not_a_knot - diag[n - 1] = dx[n - 3]; - subdiag[n - 2] = x[n - 1] - x[n - 3]; - s[n - 1] = (dd[n - 2] * dx[n - 3] * (dx[n - 2] + 2 * (x[n - 1] - x[n - 3])) - + dd[n - 3] * dx[n - 2] * dx[n - 2]) - / (x[n - 1] - x[n - 3]); + d[n - 1] = dx[n - 3]; + l[n - 2] = dx[n - 3] + dx[n - 2]; + dy[n - 1] = (dd[n - 2] * dx[n - 3] * (dx[n - 2] + 2 * l[n - 2]) + + dd[n - 3] * dx[n - 2] * dx[n - 2]) / l[n - 2]; } - int NRHS = 1; - int LDB = n; - int INFO = 0; - int N = n; - - dgtsv_(&N, &NRHS, subdiag.data(), diag.data(), supdiag.data(), s, &LDB, &INFO); + int nrhs = 1; + int ldb = n; + int info = 0; + dgtsv_(&n, &nrhs, l, d, u, dy, &ldb, &info); } } } -void CubicSpline::build(const int n, - const double* const x, - const double* const y, - BoundaryCondition bc_start, - BoundaryCondition bc_end, - const double deriv_start, - const double deriv_end) -{ - cleanup(); - - n_ = n; - x_ = new double[n]; - y_ = new double[n]; - s_ = new double[n]; - - std::memcpy(x_, x, sizeof(double) * n); - std::memcpy(y_, y, sizeof(double) * n); - build(n_, x_, y_, s_, bc_start, bc_end, deriv_start, deriv_end); - - double dx_avg = (x[n - 1] - x[0]) / (n - 1); - is_uniform_ = std::all_of(x, x + n, - [dx_avg, &x](const double& xi) -> bool { return std::abs(xi - (&xi - x) * dx_avg) < 1e-15; }); -} - -void CubicSpline::eval(const int n, - const double* const x, - const double* const y, - const double* const s, - const int n_interp, - const double* const x_interp, - double* const y_interp, - double* const dy_interp) -{ - check_interp(n, x, y, s, n_interp, x_interp, y_interp, dy_interp); - _eval(x, y, s, n_interp, x_interp, y_interp, dy_interp, _gen_search(n, x)); -} - -void CubicSpline::eval(const int n_interp, - const double* const x_interp, - double* const y_interp, - double* const dy_interp) -{ - check_interp(n_, x_, y_, s_, n_interp, x_interp, y_interp, dy_interp); - _eval(x_, y_, s_, n_interp, x_interp, y_interp, dy_interp, _gen_search(n_, x_, is_uniform_)); -} -std::function CubicSpline::_gen_search(const int n, const double* const x, int is_uniform) +void CubicSpline::_cubic( + int n, + const double* w, + const double* c0, + const double* c1, + const double* c2, + const double* c3, + double* y, + double* dy, + double* d2y +) { - if (is_uniform != 0 && is_uniform != 1) + if (y) { - double dx_avg = (x[n - 1] - x[0]) / (n - 1); - is_uniform = std::all_of(x, x + n, - [dx_avg, &x] (const double& xi) { return std::abs(xi - (&xi - x) * dx_avg) < 1e-15; }); + for (int i = 0; i < n; ++i) + { + y[i] = ((c3[i] * w[i] + c2[i]) * w[i] + c1[i]) * w[i] + c0[i]; + } } - if (is_uniform) + if (dy) { - double dx = x[1] - x[0]; - return [dx, n, x](double xi) -> int { return xi == x[n - 1] ? n - 2 : xi / dx; }; + for (int i = 0; i < n; ++i) + { + dy[i] = (3.0 * c3[i] * w[i] + 2.0 * c2[i]) * w[i] + c1[i]; + } } - else + + if (d2y) { - return [n, x](double xi) -> int { return (std::upper_bound(x, x + n, xi) - x) - 1; }; + for (int i = 0; i < n; ++i) + { + d2y[i] = 6.0 * c3[i] * w[i] + 2.0 * c2[i]; + } } } -void CubicSpline::_eval(const double* const x, - const double* const y, - const double* const s, - const int n_interp, - const double* const x_interp, - double* const y_interp, - double* const dy_interp, - std::function search) -{ - void (*poly_eval)(double, double, double, double, double, double*, double*) = - y_interp ? (dy_interp ? _poly_eval: _poly_eval) : _poly_eval; - - for (int i = 0; i != n_interp; ++i) - { - int p = search(x_interp[i]); - double w = x_interp[i] - x[p]; - - double dx = x[p + 1] - x[p]; - double dd = (y[p + 1] - y[p]) / dx; - - double c0 = y[p]; - double c1 = s[p]; - double c3 = (s[p] + s[p + 1] - 2.0 * dd) / (dx * dx); - double c2 = (dd - s[p]) / dx - c3 * dx; - poly_eval(w, c0, c1, c2, c3, y_interp + i, dy_interp + i); - } +int CubicSpline::_index(int n, const double* knots, double x) +{ + int i = (std::upper_bound(knots, knots + n, x) - knots) - 1; + return i - (i == n - 1); } -void CubicSpline::solve_cyctri(const int n, - double* const diag, - double* const supdiag, - double* const subdiag, - double* const b) + +int CubicSpline::_index(int n, double x0, double dx, double x) { - // This function uses the Sherman-Morrison formula to convert a cyclic tridiagonal linear system - // into a tridiagonal one. + int i = (x - x0) / dx; + return i - (i == n - 1); +} - // NOTE all diags will be overwritten in this function! - // flexible non-zero parameters that can affect the condition number of the tridiagonal linear system below - // may have some smart choice, set to 1 for now +void CubicSpline::_solve_cyctri(int n, double* d, double* u, double* l, double* b) +{ + // flexible non-zero parameters that can affect the condition number of the + // tridiagonal linear system double alpha = 1.0; - double beta = 1.0; + double beta = -d[0] / u[n - 1]; - double* bp = new double[2 * n]; - std::memcpy(bp, b, sizeof(double) * n); + std::vector bp(2 * n, 0.0); + std::copy(b, b + n, bp.begin()); bp[n] = 1. / alpha; bp[2 * n - 1] = 1. / beta; - for (int i = n + 1; i != 2 * n - 1; ++i) - { - bp[i] = 0.0; - } - diag[0] -= supdiag[n - 1] * beta / alpha; - diag[n - 1] -= subdiag[n - 1] * alpha / beta; + d[0] -= u[n - 1] * beta / alpha; + d[n - 1] -= l[n - 1] * alpha / beta; int nrhs = 2; int info = 0; - int N = n; int ldb = n; - dgtsv_(&N, &nrhs, subdiag, diag, supdiag, bp, &ldb, &info); + dgtsv_(&n, &nrhs, l, d, u, bp.data(), &ldb, &info); - double fac = (beta * supdiag[n - 1] * bp[0] + alpha * subdiag[n - 1] * bp[n - 1]) - / (1. + beta * supdiag[n - 1] * bp[n] + alpha * subdiag[n - 1] * bp[2 * n - 1]); + double fac = (beta * u[n - 1] * bp[0] + alpha * l[n - 1] * bp[n - 1]) + / (1. + beta * u[n - 1] * bp[n] + alpha * l[n - 1] * bp[2 * n - 1]); - std::memcpy(b, bp, sizeof(double) * n); - for (int i = 0; i != n; ++i) - { - b[i] -= fac * bp[n + i]; - } - - delete[] bp; + std::transform(bp.begin(), bp.begin() + n, bp.begin() + n, b, + [fac](double yi, double zi) { return yi - fac * zi; }); } -} // namespace ModuleBase + diff --git a/source/module_base/cubic_spline.h b/source/module_base/cubic_spline.h index c025dd7b6e..92936704de 100644 --- a/source/module_base/cubic_spline.h +++ b/source/module_base/cubic_spline.h @@ -1,82 +1,136 @@ -#ifndef CUBIC_SPLINE_INTERPOLATOR_H_ -#define CUBIC_SPLINE_INTERPOLATOR_H_ +#ifndef CUBIC_SPLINE_INTERPOLATION_H_ +#define CUBIC_SPLINE_INTERPOLATION_H_ -#include "lapack_connector.h" +#include +#include namespace ModuleBase { -/*! - * @brief A class that performs cubic spline interplations. +/** + * @brief Cubic spline interplation. * - * This class interpolates a given set of data points (x[i],y[i]) (i=0,...,n-1) - * by piecewise cubic polynomials with continuous first and second derivatives - * at x[i]. + * Interpolating a set of data points (x[i], y[i]) (i=0,...,n-1) by piecewise + * cubic polynomials * - * There are two ways to use this class. The first way treats the class as an - * interpolant object, and the second way uses the static member functions. + * p_i(x) = c0[i] + c1[i]*(x-x[i]) + c2[i]*(x-x[i])^2 + c3[i]*(x-x[i])^3 * - * Usage-1: interpolant object + * with continuous first and second derivatives. (p_i(x) is defined on the + * interval [x[i], x[i+1]]) * - * CubicSpline cubspl; + * There are two ways to use this class. The first way treats class objects as + * interpolants; the second way uses static member functions. * - * // build the interpolant - * // n is the number of data points (x[i],y[i]) (i=0,...,n-1) - * // by default "not-a-knot" boundary condition is used at both ends - * cubspl.build(n, x, y); + * Usage-1: object as interpolant * - * // alternatively one can specify the first or second derivatives at the ends - * cubspl.build(n, x, y, CubicSpline::BoundaryCondition::second_deriv, - * CubicSpline::BoundaryCondition::second_deriv, 1.0, -2.0); + * //--------------------------------------------------------------------- + * // basic usage + * //--------------------------------------------------------------------- + * // build the interpolant object + * // n is the number of data points (x[i], y[i]) (i=0,...,n-1) + * CubicSpline cubspl(n, x, y); * - * // not-a-knot, first_deriv & second_deriv can be independently applied to - * // each end; two ends do not necessarily have the same type of condition - * cubspl.build(n, x, y, CubicSpline::BoundaryCondition::first_deriv, - * CubicSpline::BoundaryCondition::not-a-knot, 1.0); + * // evaluates the interpolant at multiple places (x_interp) + * // n_interp is the number of places to evaluate the interpolant + * cubspl.eval(n_interp, x_interp, y_interp); // values are returned in y_interp * - * // periodic boundary condition needs to be specified at both ends - * // and y[0] must equal y[n-1] in this case - * cubspl.build(n, x, y, CubicSpline::BoundaryCondition::periodic, - * CubicSpline::BoundaryCondition::periodic ); + * // evaluates both the values and first derivatives at x_interp + * cubspl.eval(n_interp, x_interp, y_interp, dy_interp); * - * // calculate the values of interpolant at x_interp[i] - * cubspl.eval(n_interp, x_interp, y_interp); + * // evaluates the second derivative only + * cubspl.eval(n_interp, x_interp, nullptr, nullptr, d2y_interp); + * + * //--------------------------------------------------------------------- + * // evenly spaced knots + * //--------------------------------------------------------------------- + * // Interpolants with evenly spaced knots can be built by a different + * // constructor, which allows faster evaluation due to quicker index lookup. + * + * // build an interpolant with evenly spaced knots x[i] = x0 + i*dx + * CubicSpline cubspl(n, x0, dx, y); + * + * //--------------------------------------------------------------------- + * // boundary conditions + * //--------------------------------------------------------------------- + * // By default "not-a-knot" boundary condition is applied to both ends. + * // Other supported boundary conditions include first/second derivatives + * // and periodic boundary condition. + * + * // build an interpolant with f''(start) = 1.0 and f'(end) = 3.0 + * CubicSpline cubspl(n, x, y, + * {CubicSpline::BoundaryType::second_deriv, 1.0}, + * {CubicSpline::BoundaryType::first_deriv, 3.0}); + * + * // build an interpolant with periodic boundary condition + * CubicSpline cubspl(n, x, y, // y[0] must equal y[n-1] + * {CubicSpline::BoundaryType::periodic}, + * {CubicSpline::BoundaryType::periodic}); + * + * //--------------------------------------------------------------------- + * // multiple interpolants + * //--------------------------------------------------------------------- + * // Once an object is constructed, more interpolants that share the same + * // knots can be added. + * // Such interpolants can be evaluated simultaneously at a single place. + * + * // build an object with 5 interpolants + * CubicSpline cubspl5(n, x, y); + * cubspl5.reserve(5); // reduce memory reallocations & data copies + * cubspl5.add(y2); + * cubspl5.add(y3, {CubicSpline::BoundaryType::first_deriv, 1.0}, {}); + * cubspl5.add(y4, {}, {CubicSpline::BoundaryType::second_deriv, 2.0}); + * cubspl5.add(y5); + * + * // evaluates the five interpolants simultaneously at a single place + * cubspl5.multi_eval(x_interp, y_interp) + * + * // evaluate the first and third interpolants at a single place + * std::vector ind = {0, 2}; + * cubspl5.multi_eval(ind.size(), ind.data(), x_interp, y_interp) + * + * // evaluate the last interpolant (i_spline == 4) at multiple places + * cubspl5.eval(n_interp, x_interp, y_interp, nullptr, nullptr, 4) * - * // calculate the values & derivatives of interpolant at x_interp[i] - * cubspl.eval(n_interp, x_interp, y_interp, dy_interp); * * Usage-2: static member functions * - * // gets the first-derivatives (s) at knots - * // may build with various boundary conditions as above - * CubicSpline::build(n, x, y, s); + * // step-1: computes the first-derivatives at knots + * // boundary conditions defaulted to "not-a-knot" + * CubicSpline::build(n, x, y, {}, {}, dy); * - * // evaluates the interpolant with knots, values & derivatives - * CubicSpline::eval(n, x, y, s, n_interp, x_interp, y_interp, dy_interp); + * // Various boundary conditions and evenly spaced knots are supported + * // in the same way as the interpolant object. * - * */ + * // step-2: computes the interpolated values & derivatives + * CubicSpline::eval(n, x, y, dy, n_interp, x_interp, y_interp, dy_interp); + * + * // Simultaneous evaluation of multiple interpolants are not supported + * // for static functions. + * + */ class CubicSpline -{ - public: - CubicSpline(){}; - CubicSpline(CubicSpline const&); - CubicSpline& operator=(CubicSpline const&); - - ~CubicSpline() { cleanup(); } - - /*! - * @brief Boundary conditions of cubic spline interpolations. - * - * Available boundary conditions include: - * - not_a_knot: the first two pieces at the start or the last two at the end - * are the same polynomial, i.e., x[1] or x[n-2] is not a "knot"; - * - first_deriv: User-defined first derivative; - * - second_deriv: User-defined second derivative; - * - periodic: The first and second derivatives at two ends are continuous. - * This condition requires that y[0] = y[n-1] and it must be - * applied to both ends. - * */ - enum class BoundaryCondition +{ + //***************************************************************** + // boundary condition + //***************************************************************** + +public: + + /** + * @brief Types of cubic spline boundary conditions. + * + * Supported types include: + * - not_a_knot The first or last two pieces are the same polynomial, + * i.e., x[1] or x[n-2] is not a "knot". This does not + * rely on any prior knowledge of the original function + * and is the default option. + * - first_deriv user-defined first derivative + * - second_deriv user-defined second derivative + * - periodic the first and second derivatives at two ends are continuous. + * This condition requires that y[0] = y[n-1] and it must be + * applied to both ends + */ + enum class BoundaryType { not_a_knot, first_deriv, @@ -84,189 +138,436 @@ class CubicSpline periodic }; - /*! - * @brief Builds the interpolant. - * - * By calling this function, the class object computes and stores the first-order - * derivatives at knots by from given data points and boundary conditions, which - * makes the object an interpolant. - * */ - void build(const int n, //!< [in] number of data points - const double* const x, //!< [in] x coordiates of data points, must be strictly increasing - const double* const y, //!< [in] y coordiates of data points - BoundaryCondition bc_start = BoundaryCondition::not_a_knot, //!< [in] boundary condition at start - BoundaryCondition bc_end = BoundaryCondition::not_a_knot, //!< [in] boundary condition at end - const double deriv_start = 0.0, //!< [in] first or second derivative at the start, - //!< used if bc_start is first_deriv or second_deriv - const double deriv_end = 0.0 //!< [in] first or second derivative at the end, - //!< used if bc_end is first_deriv or second_deriv - ); - /*! - * @brief Evaluates the interpolant. + /** + * @brief Boundary condition for cubic spline interpolation. * - * This function evaluates the interpolant at x_interp[i]. - * On finish, interpolated values are placed in y_interp, - * and the derivatives at x_interp[i] are placed in dy_interp. + * An object of this struct represents an actual boundary condition at one end, + * which contains a type and possibly a value (first/second_deriv only). * - * If y_interp or dy_interp is nullptr, the corresponding values are - * not calculated. They must not be nullptr at the same time. + */ + struct BoundaryCondition + { + // for not_a_knot and periodic + BoundaryCondition(BoundaryType type = BoundaryType::not_a_knot); + + // for first/second_deriv + BoundaryCondition(BoundaryType type, double val); + + BoundaryType type; + double val = 0.0; + }; + + + //***************************************************************** + // interpolant object + //***************************************************************** + +public: + + CubicSpline() = delete; + CubicSpline(CubicSpline const&) = default; + CubicSpline(CubicSpline &&) = default; + + CubicSpline& operator=(CubicSpline const&) = default; + CubicSpline& operator=(CubicSpline &&) = default; + + ~CubicSpline() = default; + + + /** + * @brief Builds an interpolant object. * - * @note the interpolant must be built before calling this function. - * */ - void eval(const int n, //!< [in] number of points to evaluate the interpolant - const double* const x_interp, //!< [in] places where the interpolant is evaluated; - //!< must be within [x_[0], x_[n-1]] - double* const y_interp, //!< [out] interpolated values - double* const dy_interp = nullptr //!< [out] derivatives at x_interp + * Constructing a cubic spline interpolant from a set of data points + * (x[i], y[i]) (i=0,1,...,n-1) and boundary conditions. + * + * @param[in] n number of data points + * @param[in] x x coordinates of data points + * ("knots", must be strictly increasing) + * @param[in] y y coordinates of data points + * @param[in] bc_start boundary condition at start + * @param[in] bc_end boundary condition at end + * + */ + CubicSpline( + int n, + const double* x, + const double* y, + const BoundaryCondition& bc_start = {}, + const BoundaryCondition& bc_end = {} ); - /// knots of the interpolant - const double* x() const { return x_; } - - /// values at knots - const double* y() const { return y_; } - - /// first-order derivatives at knots - const double* s() const { return s_; } - - static void build(const int n, //!< [in] number of data points - const double* const x, //!< [in] x coordiates of data points, must be strictly increasing - const double* const y, //!< [in] y coordiates of data points - double* const s, //!< [out] first-order derivatives at knots - BoundaryCondition bc_start = BoundaryCondition::not_a_knot, //!< [in] boundary condition at start - BoundaryCondition bc_end = BoundaryCondition::not_a_knot, //!< [in] boundary condition at end - const double deriv_start = 0.0, //!< [in] first or second derivative at the start, - //!< used if bc_start is first_deriv or second_deriv - const double deriv_end = 0.0 //!< [in] first or second derivative at the end, - //!< used if bc_end is first_deriv or second_deriv + + /** + * @brief Builds an interpolant object with evenly-spaced knots. + * + * Constructing a cubic spline interpolant from a set of data points + * (x0+i*dx, y[i]) (i=0,1,...,n-1) and boundary conditions. + * + * @param[in] n number of data points + * @param[in] x0 x coordinate of the first data point (first knot) + * @param[in] dx spacing between knots (must be positive) + * @param[in] y y coordinates of data points + * @param[in] bc_start boundary condition at start + * @param[in] bc_end boundary condition at end + * + */ + CubicSpline( + int n, + double x0, + double dx, + const double* y, + const BoundaryCondition& bc_start = {}, + const BoundaryCondition& bc_end = {} ); - static void eval(const int n, //!< [in] number of knots - const double* const x, //!< [in] knots of the interpolant - const double* const y, //!< [in] values at knots - const double* const s, //!< [in] first-order derivatives at knots - const int n_interp, //!< [in] number of points to evaluate the interpolant - const double* const x_interp, //!< [in] places where the interpolant is evaluated; - //!< must be within [x_[0], x_[n-1]] - double* const y_interp, //!< [out] interpolated values - double* const dy_interp = nullptr //!< [out] derivatives at x_interp + + /** + * @brief Add an interpolant that shares the same knots. + * + * An object of this class can hold multiple interpolants with the same knots. + * Once constructed, more interpolants sharing the same knots can be added by + * this function. Multiple interpolants can be evaluated simultaneously at a + * single place by multi_eval. + * + * @param[in] y y coordinates of data points + * @param[in] bc_start boundary condition at start + * @param[in] bc_end boundary condition at end + * + */ + void add( + const double* y, + const BoundaryCondition& bc_start = {}, + const BoundaryCondition& bc_end = {} ); - private: - //! number of data points + + /** + * @brief Evaluates a single interpolant at multiple places. + * + * @param[in] n_interp number of places to evaluate the interpolant + * @param[in] x_interp places where an interpolant is evaluated + * (must be within the range of knots) + * @param[out] y_interp values at x_interp + * @param[out] dy_interp first derivatives at x_interp + * @param[out] d2y_interp second derivatives at x_interp + * @param[in] i_spline index of the interpolant to evaluate + * + * @note pass nullptr to any of the output would suppress the corresponding calculation + * + */ + void eval( + int n_interp, + const double* x_interp, + double* y_interp, + double* dy_interp = nullptr, + double* d2y_interp = nullptr, + int i_spline = 0 + ) const; + + + /** + * @brief Evaluates multiple interpolants at a single place. + * + * @param[in] n_spline number of interpolants to evaluate + * @param[in] i_spline indices of interpolants to evaluate + * @param[in] x_interp place where interpolants are evaluated + * (must be within the range of knots) + * @param[out] y_interp values at x_interp + * @param[out] dy_interp first derivatives at x_interp + * @param[out] d2y_interp second derivatives at x_interp + * + * @note pass nullptr to any of the output would suppress the corresponding calculation + * + */ + void multi_eval( + int n_spline, + const int* i_spline, + double x_interp, + double* y_interp, + double* dy_interp = nullptr, + double* d2y_interp = nullptr + ) const; + + + /** + * @brief Evaluates all interpolants at a single place. + * + * @param[in] x_interp place where interpolants are evaluated + * (must be within the range of knots) + * @param[out] y_interp values at x_interp + * @param[out] dy_interp first derivatives at x_interp + * @param[out] d2y_interp second derivatives at x_interp + * + * @note pass nullptr to any of the output would suppress the corresponding calculation + * + */ + void multi_eval( + double x_interp, + double* y_interp, + double* dy_interp = nullptr, + double* d2y_interp = nullptr + ) const; + + + /** + * @brief Reserves memory for holding more interpolants. + * + * By default this class does not reserve memory for multiple interpolants. + * Without reservation, whenever a new interpolant is added, memory has to + * be reallocated and old data copied, which could waste a lot of time if + * there's a large number of interpolants to add. + * + * This function help avoid repetitive memory reallocations and data copies + * by a one-shot reservation. + * + * @param[in] n_spline expected total number of interpolants + * + */ + void reserve(int n_spline) { y_.reserve(n_spline * n_ * 2); } + + + /// heap memory usage in bytes + size_t heap_usage() const { return (x_.capacity() + y_.capacity()) * sizeof(double); } + + /// first knot + double xmin() const { return xmin_; } + + /// last knot + double xmax() const { return xmax_; } + + +private: + + /// number of cubic spline interpolants + int n_spline_ = 0; + + /// number of knots int n_ = 0; - //! knots (x coordinates of data points) - double* x_ = nullptr; - - //! values at knots - double* y_ = nullptr; - - //! first-order derivatives at knots - double* s_ = nullptr; - - //! A flag that tells whether the knots are evenly spaced. - bool is_uniform_ = false; - - //! Numerical threshold for determining whether the knots are evenly spaced. - /*! - * The knots are considered uniform (evenly spaced) if for every i from 0 to n-2 - * - * abs( (x[i+1]-x[i]) - (x[n-1]-x[0])/(n-1) ) < uniform_thr_ - * */ - double uniform_thr_ = 1e-15; - - /// Checks whether the input arguments are valid for building a cubic spline. - static void check_build(const int n, - const double* const x, - const double* const y, - BoundaryCondition bc_start, - BoundaryCondition bc_end); - - /// Checks whether the input arguments are valid for evaluating a cubic spline. - static void check_interp(const int n, - const double* const x, - const double* const y, - const double* const s, - const int n_interp, - const double* const x_interp, - double* const y_interp, - double* const dy_interp); - - //! Solves a cyclic tridiagonal linear system. - /*! - * This function solves a cyclic tridiagonal linear system A*x=b where b - * is a vector and A is given by - * - * D[0] U[0] L[n-1] - * L[0] D[1] U[1] - * L[1] D[2] U[2] - * ... ... ... - * L[n-3] D[n-2] U[n-2] - * U[n-1] L[n-2] D[n-1] - * - * On finish, b is overwritten by the solution. - * - * Sherman-Morrison formula is used to convert the problem into a tridiagonal - * linear system, after which the problem can be solved by dgtsv efficiently. - * - * @note D, L, U are all modified in this function, so use with care! - * */ - static void solve_cyctri(const int n, //!< [in] size of the linear system - double* const D, //!< [in] main diagonal - double* const U, //!< [in] superdiagonal - double* const L, //!< [in] subdiagonal - double* const b //!< [in,out] right hand side of the linear system; - //!< will be overwritten by the solution on finish. + /// first knot + double xmin_ = 0.0; + + /// last knot + double xmax_ = 0.0; + + /// spacing between knots (only used for evenly-spaced knots) + double dx_ = 0.0; + + /// knots of the spline polynomial (remains empty for evenly-spaced knots) + std::vector x_; + + /// values and first derivatives at knots + std::vector y_; + + + //***************************************************************** + // static functions + //***************************************************************** + +public: + + /** + * @brief Computes the first derivatives at knots for cubic spline + * interpolation. + * + * @param[in] n number of data points + * @param[in] x x coordinates of data points + * ("knots", must be strictly increasing) + * @param[in] y y coordinates of data points + * @param[in] bc_start boundary condition at start + * @param[in] bc_end boundary condition at end + * @param[out] dy first derivatives at knots + * + */ + static void build( + int n, + const double* x, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end, + double* dy ); - //! Wipes off the interpolant (if any) and deallocates memories. - void cleanup(); - /// Evaluates a cubic polynomial and its derivative. - template - static void _poly_eval(double w, double c0, double c1, double c2, double c3, double* y, double* dy) - { - if (EvalY) - { - *y = ((c3 * w + c2) * w + c1) * w + c0; - } - - if (EvalDy) - { - *dy = (3.0 * c3 * w + 2.0 * c2) * w + c1; - } - } - - /*! - * @brief Generates a function that returns the index of the left knot of the - * spline polynomial to be evaluated. - * - * This function takes the knots of the interpolant and returns a function that - * takes a value x_interp and returns the index of the left knot of the spline - * polynomial. If "is_uniform" is not 0/1, this function will checks whether - * the knots are evenly spaced. The returned function makes use of "is_uniform" - * to speed up the search. - * */ - static std::function _gen_search(const int n, const double* const x, int is_uniform = -1); - - /*! - * @brief Evaluates a cubic spline with given knots, values and derivatives. - * */ - static void _eval(const double* const x, //!< [in] knots of the interpolant - const double* const y, //!< [in] values at knots - const double* const s, //!< [in] first-order derivatives at knots - const int n_interp, //!< [in] number of points to evaluate the interpolant - const double* const x_interp, //!< [in] places where the interpolant is evaluated; - //!< must be within [x_[0], x_[n-1]] - double* const y_interp, //!< [out] interpolated values - double* const dy_interp, //!< [out] derivatives at x_interp - std::function search //!< [in] a function that returns the index of the left - // knot of the spline polynomial to be evaluated + /** + * @brief Computes the first derivatives at evenly-spaced knots for + * cubic spline interpolation. + * + * @param[in] n number of data points + * @param[in] dx spacing between knots (must be positive) + * @param[in] y y coordinates of data points + * @param[in] bc_start boundary condition at start + * @param[in] bc_end boundary condition at end + * @param[out] dy first derivatives at knots + * + */ + static void build( + int n, + double dx, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end, + double* dy + ); + + + /** + * @brief Evaluates a cubic spline polynomial at multiple places. + * + * @param[in] n number of knots + * @param[in] x knots (must be strictly increasing) + * @param[in] y values at knots + * @param[in] dy first derivatives at knots + * @param[in] n_interp number of places to evaluate the interpolant + * @param[in] x_interp places where the interpolant is evaluated + * (must be within the range of knots) + * @param[out] y_interp values at x_interp + * @param[out] dy_interp first derivatives at x_interp + * @param[out] d2y_interp second derivatives at x_interp + * + * @note pass nullptr to any of the output would suppress the corresponding calculation + * + */ + static void eval( + int n, + const double* x, + const double* y, + const double* dy, + int n_interp, + const double* x_interp, + double* y_interp, + double* dy_interp = nullptr, + double* d2y_interp = nullptr + ); + + + /** + * @brief Evaluates a cubic spline polynomial with evenly spaced knots. + * + * @param[in] n number of knots + * @param[in] x0 first knot + * @param[in] dx spacing between knots + * @param[in] y values at knots + * @param[in] dy first derivatives at knots + * @param[in] n_interp number of places to evaluate the interpolant + * @param[in] x_interp places where the interpolant is evaluated + * (must be within the range of knots) + * @param[out] y_interp values at x_interp + * @param[out] dy_interp first derivatives at x_interp + * @param[out] d2y_interp second derivatives at x_interp + * + * @note pass nullptr to any of the output would suppress the corresponding calculation + * + */ + static void eval( + int n, + double x0, + double dx, + const double* y, + const double* dy, + int n_interp, + const double* x_interp, + double* y_interp, + double* dy_interp = nullptr, + double* d2y_interp = nullptr + ); + + +private: + + /// Computational routine for building cubic spline interpolant + static void _build( + int n, + const double* dx, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end, + double* dy + ); + + + /** + * @brief Segment index lookup. + * + * Given a strictly increasing array x and a target within the range of + * x, this function returns an index i such that x[i] <= target < x[i+1] + * if target != x[n-1], or n-2 if t == x[n-1]. + * + */ + static inline int _index(int n, const double* x, double target); + + + /// Segment index lookup (evenly spaced knots). + static inline int _index(int n, double x0, double dx, double target); + + + /// Evaluates a batch of cubic polynomials. + static inline void _cubic( + int n, + const double* w, + const double* c0, + const double* c1, + const double* c2, + const double* c3, + double* y, + double* dy, + double* d2y + ); + + + /// Asserts that the input arguments are valid for constructing a cubic spline. + static void _validate_build( + int n, + const double* dx, + const double* y, + const BoundaryCondition& bc_start, + const BoundaryCondition& bc_end ); + + + /// Asserts that the input arguments are valid for interpolating a cubic spline. + static void _validate_eval( + int n, + const double (&u)[2], + const double* x, + const double* y, + const double* dy, + int n_interp, + const double* x_interp + ); + + + /** + * @brief Solves a cyclic tridiagonal linear system. + * + * A cyclic tridiagonal linear system A*x=b where b is a vector and + * + * -- -- + * | d[0] u[0] l[n-1] | + * | l[0] d[1] u[1] | + * A = | l[1] d[2] u[2] | + * | ... ... ... | + * | l[n-3] d[n-2] u[n-2] | + * | u[n-1] l[n-2] d[n-1] | + * -- -- + * + * is transformed to a tridiagonal linear system by the Sherman-Morrison + * formula, and then solved by dgtsv. + * + * @param[in] n size of the linear system + * @param[in] d main diagonal + * @param[in] u superdiagonal + * @param[in] l subdiagonal + * @param[in,out] b right hand side of the linear system; will be + * overwritten by the solution on finish. + * + * @note d, l, u are all overwritten in this function. + * + */ + static void _solve_cyctri(int n, double* d, double* u, double* l, double* b); }; -}; // namespace ModuleBase +} // namespace ModuleBase #endif diff --git a/source/module_base/test/cubic_spline_test.cpp b/source/module_base/test/cubic_spline_test.cpp index 9c15d98230..9eb4e8f05d 100644 --- a/source/module_base/test/cubic_spline_test.cpp +++ b/source/module_base/test/cubic_spline_test.cpp @@ -1,1468 +1,459 @@ #include "module_base/cubic_spline.h" #include +#include +#include +#include +#include #include "gtest/gtest.h" -#include "module_base/constants.h" -#include "stdio.h" - -#ifdef __MPI -#include -#endif using ModuleBase::CubicSpline; -using ModuleBase::PI; - -/*********************************************************** - * Unit test of class CubicSpline - ***********************************************************/ +using BoundaryCondition = CubicSpline::BoundaryCondition; +using BoundaryType = CubicSpline::BoundaryType; -/*! Tested functions: +/** + * @brief Unit test of class CubicSpline + * + * Tested functions include: * * - build - * Constructs the cubic spline interpolant from given - * data points and boundary condition. + * Constructs a cubic spline interpolant from + * a set of data points and boundary conditions. * * - eval - * Evaluates the interpolant at certain points. - * */ -class CubicSplineTest : public ::testing::Test -{ - protected: - /// Allocates buffers - void SetUp(); - - /// Deallocates buffers - void TearDown(); - - int n_max = 20; ///< size of each buffer - double* x_; ///< buffer for the x coordinates of data points - double* y_; ///< buffer for the y coordinates of data points - double* s_; ///< buffer for the first-derivatives of the interpolant - double* x_interp_; ///< places where the interpolant is evaluated - double* y_interp_; ///< values of the interpolant at x_interp_ - double* y_ref_; ///< reference values for y_interp_ drawn from - ///< scipy.interpolate.CubicSpline - - double tol = 1e-15; ///< tolerance for element-wise numerical error -}; - -void CubicSplineTest::SetUp() -{ - x_ = new double[n_max]; - y_ = new double[n_max]; - s_ = new double[n_max]; - x_interp_ = new double[n_max]; - y_interp_ = new double[n_max]; - y_ref_ = new double[n_max]; -} - -void CubicSplineTest::TearDown() -{ - delete[] x_; - delete[] y_; - delete[] s_; - delete[] x_interp_; - delete[] y_interp_; - delete[] y_ref_; -} -#define MAX(x, y) ((x) > (y) ? (x) : (y)) - -/** - * @brief + * Evaluates a single interpolant at multiple places. * - * @param x:theInterpolation point group - * @param func:Fourth derivative of the interpolating point function - * @param n:Interpolating point set length - * @return double:Theoretical estimation error - */ -double count_err(double* x, std::function func, int n) -{ - double maxf4x = -10000000000; - double maxh = 0; - - for (int i = 0; i != n; ++i) - { - maxf4x = MAX(abs(func(x[i])), maxf4x); - } - // printf("maxf4x = %.8lf\n", maxf4x); - for (int i = 1; i != n; ++i) - { - maxh = MAX(x[i + 1] - x[i], maxh); - } - - // if(func==2) printf("maxh = %.8lf,maxf4x = %.8lf\n", maxh,maxf4x); - double err = (5.0 / 384.0) * maxf4x * pow(maxh, 4); - return MAX(err, 1e-15); - ; -} - -TEST_F(CubicSplineTest, NotAKnot) -{ - CubicSpline cubspl; - - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = std::sqrt(i); - y_[i] = std::sin(x_[i]); - } - - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.57; - x_interp_[3] = 2.00; - x_interp_[4] = 2.99; - x_interp_[5] = 3.00; - - y_ref_[0] = 0.; - y_ref_[1] = 0.0105903444284005; - y_ref_[2] = 1.0000633795463434; - y_ref_[3] = 0.9092974268256817; - y_ref_[4] = 0.1510153464180796; - y_ref_[5] = 0.1411200080598672; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } - - // static member function - ModuleBase::CubicSpline::build(n, - x_, - y_, - s_, - ModuleBase::CubicSpline::BoundaryCondition::not_a_knot, - ModuleBase::CubicSpline::BoundaryCondition::not_a_knot); - ModuleBase::CubicSpline::eval(n, x_, y_, s_, ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } -} - -TEST_F(CubicSplineTest, PeriodicAndUniform) -{ - CubicSpline cubspl; - - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = i * 2 * PI / (n - 1); - y_[i] = std::cos(x_[i]); - } - - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::periodic, CubicSpline::BoundaryCondition::periodic); - - int ni = 5; - x_interp_[0] = 0.0; - x_interp_[1] = 0.5 * PI; - x_interp_[2] = PI; - x_interp_[3] = 1.5 * PI; - x_interp_[4] = 2 * PI; - - y_ref_[0] = 1.0000000000000000e+00; - y_ref_[1] = 1.4356324132368183e-04; - y_ref_[2] = -9.9930291851807085e-01; - y_ref_[3] = 1.4356324132349871e-04; - y_ref_[4] = 1.0000000000000000e+00; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } -} - -TEST_F(CubicSplineTest, FirstDeriv) -{ - CubicSpline cubspl; - - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = std::sqrt(i); - y_[i] = std::exp(-x_[i]); - } - - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - -3, - -0.5); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 2.0; - x_interp_[4] = 2.54; - x_interp_[5] = 3.00; - - y_ref_[0] = 1.; - y_ref_[1] = 0.9704131180863818; - y_ref_[2] = 0.1367376505691157; - y_ref_[3] = 0.1353352832366127; - y_ref_[4] = 0.0798871927951471; - y_ref_[5] = 0.0497870683678639; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } -} - -TEST_F(CubicSplineTest, TwoPoints) -{ - CubicSpline cubspl; - - int ni = 5; - x_interp_[0] = 0.0; - x_interp_[1] = 0.1; - x_interp_[2] = 0.33; - x_interp_[3] = 0.9; - x_interp_[4] = 1.0; - - x_[0] = 0; - x_[1] = 1; - y_[0] = 2.33; - y_[1] = 4.33; - - // first derivative - cubspl.build(2, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - 0.8, - 1.5); - - y_ref_[0] = 2.33; - y_ref_[1] = 2.4373; - y_ref_[2] = 2.8487171; - y_ref_[3] = 4.159700000000001; - y_ref_[4] = 4.33; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } - - // second derivative - cubspl.build(2, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - 0.8, - 1.5); - - y_ref_[0] = 2.33; - y_ref_[1] = 2.48245; - y_ref_[2] = 2.86725265; - y_ref_[3] = 4.074050000000001; - y_ref_[4] = 4.33; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } - - // periodic - y_[1] = y_[0]; - cubspl.build(2, x_, y_, CubicSpline::BoundaryCondition::periodic, CubicSpline::BoundaryCondition::periodic); - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], 2.33, tol); - } - - // "not-a-knot" is invalid for n=2 -} - -TEST_F(CubicSplineTest, ThreePoints) -{ - CubicSpline cubspl; - - int ni = 5; - x_interp_[0] = 0.0; - x_interp_[1] = 0.1; - x_interp_[2] = 0.33; - x_interp_[3] = 0.9; - x_interp_[4] = 1.0; - - // not-a-knot - x_[0] = 0; - x_[1] = 0.4; - x_[2] = 1.0; - - y_[0] = 1.2; - y_[1] = 2.5; - y_[2] = 4.8; - - cubspl.build(3, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - - y_ref_[0] = 1.2; - y_ref_[1] = 1.5075; - y_ref_[2] = 2.259025; - y_ref_[3] = 4.3875; - y_ref_[4] = 4.8; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } - - // periodic - y_[2] = y_[0]; - cubspl.build(3, x_, y_, CubicSpline::BoundaryCondition::periodic, CubicSpline::BoundaryCondition::periodic); - - y_ref_[0] = 1.2; - y_ref_[1] = 1.44375; - y_ref_[2] = 2.35383125; - y_ref_[3] = 1.2361111111111112; - y_ref_[4] = 1.2; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } -} - -/** - * @brief£º The following is the theoretical accuracy test code flow for cubic spline interpolation functions. + * - add + * Adds an interpolant that shares the same knots. * - * Boundary conditions for the cubic spline interpolation known cases of theoretical error follow the following - * formula: err<=5/384*(max(|f''''(x)|))*h^4 , which h = x_ (i + 1) - x_i + * - multi_eval + * Evaluates multiple interpolants at a single place. * - * The theoretical error test in the test code is then performed in the following three steps: - * 1. Initialize the cubic spline interpolation point and the test point set of true values - * 2. The theoretical error is calculated using the above formula through the known interpolation function - * and interpolation point. - * 3. Carry out cubic spline interpolation and compare with the real value whether the theoretical error is - * satisfied. + * - reserve + * Reserves memory for multiple interpolants. + * + * - heap_usage + * Returns the heap usage of the object. + * + * - xmin, xmax + * Returns the first and last knots. * */ -TEST_F(CubicSplineTest, FirstDeriv_sinx_uniform) -{ - CubicSpline cubspl; - - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = ((double)i / 10.0) * 2 * PI; - y_[i] = sin(x_[i]); - } - - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - 1, - 0.80901699); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 5.550; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.988258957900347; - y_ref_[5] = -0.669239857276262; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, FirstDeriv_sinx_Notuniform) +class CubicSplineTest : public ::testing::Test { - CubicSpline cubspl; +protected: + + CubicSplineTest(); + ~CubicSplineTest() = default; + + /// maximum number of knots + int n_max_; + + /// buffer for a cubic spline (x_, y_ & dy_) + std::vector spline_; + + /// knots (x-coordinates of data points) + double* x_; + + /// values at knots (y-coordinates of data points) + double* y_; + + /// derivatives at knots (computed when building a cubic spline) + double* dy_; + + /// maximum number of places to evaluate an interpolant + int n_interp_max_; + + /// buffer for interpolant evaluation + std::vector interp_; + + /// places to evaluate an interpolant + double* x_interp_; + + /// values and derivatives of the interpolant at x_interp_ + double* y_interp_; + double* dy_interp_; + double* d2y_interp_; + + /// reference values and derivatives + double* y_ref_; + double* dy_ref_; + double* d2y_ref_; + + /// y/dy/d2y tolerance for cross-check + const double tol_[3] = {1e-14, 1e-13, 1e-12}; + + + /** + * @brief Sample functions & derivatives in error bound check. + * + * @note Functions with vanishing 4-th derivative in an interval + * should in principle be interpolated exactly by a cubic spline. + * However, the presence of floating-point rounding errors would + * lead to some discrepancy between the interpolant and the original + * function. Such error is not covered by the error bound formula. + * Functions here should not include those kind of functions. + * + */ + std::vector>> f_; + + /// theoretical error bound for complete cubic spline + double error_bound( + int n, + const double* x, + const std::function& f, + int d = 0 + ) const; + + /// + void read( + const std::string& fname, + int& n, + double* x, + double* y, + BoundaryCondition& bc_start, + BoundaryCondition& bc_end, + int& n_interp, + double* x_interp, + double* y_interp, + double* dy_interp, + double* d2y_interp + ) const; +}; - int n = 10; - x_[0] = 0; - y_[0] = sin(x_[0]); - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i) * 2 * PI; - y_[9 - i] = sin(x_[9 - i]); - } - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - 1, - 1); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 6; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.756802495307928; - y_ref_[5] = -0.202090837026266; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); +CubicSplineTest::CubicSplineTest(): + n_max_(1000), + spline_(3 * n_max_), + x_(spline_.data()), + y_(x_ + n_max_), + dy_(y_ + n_max_), + n_interp_max_(1000), + interp_(7 * n_interp_max_), + x_interp_(interp_.data()), + y_interp_(x_interp_ + n_interp_max_), + dy_interp_(y_interp_ + n_interp_max_), + d2y_interp_(dy_interp_ + n_interp_max_), + y_ref_(d2y_interp_ + n_interp_max_), + dy_ref_(y_ref_ + n_interp_max_), + d2y_ref_(dy_ref_ + n_interp_max_), + f_{ + { + [](double x) { return std::sin(x); }, + [](double x) { return std::cos(x); }, + [](double x) { return -std::sin(x); }, + [](double x) { return -std::cos(x); }, + [](double x) { return std::sin(x); }, + }, + { + [](double x) { return std::exp(-x); }, + [](double x) { return -std::exp(-x); }, + [](double x) { return std::exp(-x); }, + [](double x) { return -std::exp(-x); }, + [](double x) { return std::exp(-x); }, + }, + { + [](double x) { return std::log(x); }, + [](double x) { return 1.0 / x; }, + [](double x) { return -1.0 / (x * x); }, + [](double x) { return 2.0 / (x * x * x); }, + [](double x) { return -6.0 / (x * x * x * x); }, + }, + } +{} + + +double CubicSplineTest::error_bound( + int n, + const double* x, + const std::function& d4f, + int d +) const +{ + std::vector buffer(n); + + std::adjacent_difference(x, x + n, buffer.begin()); + double max_dx = *std::max_element(buffer.begin() + 1, buffer.end()); + + auto d4f_abs = [&d4f](double x) { return std::abs(d4f(x)); }; + std::transform(x, x + n, buffer.begin(), d4f_abs); + double max_d4f = *std::max_element(buffer.begin(), buffer.end()); + + // See Carl de Boor, "A Practical Guide to Splines", Chapter V. + switch (d) + { + case 0: + return 5.0 / 384.0 * std::pow(max_dx, 4) * max_d4f; + case 1: + return 1.0 / 24.0 * std::pow(max_dx, 3) * max_d4f; + case 2: + return 3.0 / 8.0 * std::pow(max_dx, 2) * max_d4f; + default: + assert(false); // should not reach here + } +} + + +void CubicSplineTest::read( + const std::string& fname, + int& n, + double* x, + double* y, + BoundaryCondition& bc_start, + BoundaryCondition& bc_end, + int& n_interp, + double* x_interp, + double* y_interp, + double* dy_interp, + double* d2y_interp +) const +{ + std::ifstream ifs(fname); + assert(ifs.is_open()); + + std::string line, bc1, bc2; + + // read boundary conditions + std::getline(ifs, line); + std::stringstream ss(line); + ss >> bc1 >> bc2; + + auto bc_parse = [](const std::string& bc) + { + if (bc == "periodic") + { + return BoundaryCondition(BoundaryType::periodic); + } + if (bc == "not-a-knot") + { + return BoundaryCondition(BoundaryType::not_a_knot); + } + if (bc.find("first_deriv") != std::string::npos) + { + return BoundaryCondition(BoundaryType::first_deriv, + std::stod(bc.substr(12, std::string::npos))); + } + if (bc.find("second_deriv") != std::string::npos) + { + return BoundaryCondition(BoundaryType::second_deriv, + std::stod(bc.substr(13, std::string::npos))); + } + else + { + assert(false); + } + }; + + bc_start = bc_parse(bc1); + bc_end = bc_parse(bc2); + + double* data[6] = {x, y, x_interp, y_interp, dy_interp, d2y_interp}; + for (int i = 0; i < 6; ++i) + { + std::getline(ifs, line); + std::stringstream ss(line); + data[i] = std::copy(std::istream_iterator(ss), + std::istream_iterator(), data[i]); } + n = data[0] - x; + n_interp = data[2] - x_interp; } -TEST_F(CubicSplineTest, SecondDeriv_sinx_uniform) -{ - CubicSpline cubspl; - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = ((double)i / 10.0) * 2 * PI; - y_[i] = sin(x_[i]); - } - - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - 0, - 0.587785); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 5.550; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.988258957900347; - y_ref_[5] = -0.669239857276262; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, SecondDeriv_sinx_Notuniform) +TEST_F(CubicSplineTest, MultiEval) { - CubicSpline cubspl; + int n = 100; + double xmin = 0.1; + double xmax = 10; + double dx = (xmax - xmin) / (n - 1); - int n = 10; - x_[0] = 0; - y_[0] = sin(x_[0]); - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i) * 2 * PI; - y_[9 - i] = sin(x_[9 - i]); - } - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - 0, - 0); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 6; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.988258957900347; - y_ref_[5] = -0.279415498198926; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} + std::for_each(x_, x_ + n, [&](double& x) { x = xmin + (&x - x_) * dx; }); -TEST_F(CubicSplineTest, Periodic_sinx_uniform) -{ - CubicSpline cubspl; + std::transform(x_, x_ + n, y_, [this](double x) { return f_[0][0](x); }); + CubicSpline cubspl(n, xmin, dx, y_, + {BoundaryType::first_deriv, f_[0][1](xmin)}, + {BoundaryType::first_deriv, f_[0][1](xmax)}); - int n = 10; - for (int i = 0; i != n - 1; ++i) - { - x_[i] = ((double)i / 10.0) * 2 * PI; - y_[i] = sin(x_[i]); - } - x_[9] = 2 * PI; - y_[9] = sin(0); - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::periodic, CubicSpline::BoundaryCondition::periodic); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 6; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.988258957900347; - y_ref_[5] = -0.279415498198926; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) + cubspl.reserve(f_.size()); + for (size_t i = 1; i < f_.size(); ++i) { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); + std::transform(x_, x_ + n, y_, [this, i](double x) { return f_[i][0](x); }); + cubspl.add(y_, {BoundaryType::first_deriv, f_[i][1](xmin)}, + {BoundaryType::first_deriv, f_[i][1](xmax)}); } -} -TEST_F(CubicSplineTest, Periodic_sinx_Notuniform) -{ - CubicSpline cubspl; - - int n = 10; - x_[0] = 0; - y_[0] = sin(x_[0]); - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i) * 2 * PI; - y_[9 - i] = sin(x_[9 - i]); - } - x_[9] = 2 * PI; - y_[9] = sin(0); - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::periodic, CubicSpline::BoundaryCondition::periodic); - - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 6; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.988258957900347; - y_ref_[5] = -0.279415498198926; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) + std::vector> err_bound(f_.size(), std::vector(3)); + for (size_t i = 0; i < f_.size(); ++i) { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); + err_bound[i][0] = error_bound(n, x_, f_[i][4], 0); + err_bound[i][1] = error_bound(n, x_, f_[i][4], 1); + err_bound[i][2] = error_bound(n, x_, f_[i][4], 2); } -} -TEST_F(CubicSplineTest, NotAKnot_sinx_uniform) -{ - CubicSpline cubspl; + int n_interp = 1000; + double dx_interp = (xmax - xmin) / (n_interp - 1); + std::for_each(x_interp_, x_interp_ + n_interp, + [&](double& x) { x = (&x - x_interp_) * dx_interp + xmin; }); - int n = 10; - for (int i = 0; i != n; ++i) + for (int p = 0; p < n_interp; ++p) { - x_[i] = ((double)i / 10.0) * 2 * PI; - y_[i] = sin(x_[i]); - } + cubspl.multi_eval(x_interp_[p], y_interp_, dy_interp_, d2y_interp_); + double ytmp, dytmp, d2ytmp; + for (size_t i = 0; i < f_.size(); ++i) + { + EXPECT_LT(std::abs(y_interp_[i] - f_[i][0](x_interp_[p])), err_bound[i][0]); + EXPECT_LT(std::abs(dy_interp_[i] - f_[i][1](x_interp_[p])), err_bound[i][1]); + EXPECT_LT(std::abs(d2y_interp_[i] - f_[i][2](x_interp_[p])), err_bound[i][2]); - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 5; // if x=5.5, not pass - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.988258957900347; - y_ref_[5] = -0.958924274663138; - - // printf("err=%.8lf\n",err); - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); + cubspl.eval(1, &x_interp_[p], &ytmp, &dytmp, &d2ytmp, i); + EXPECT_NEAR(ytmp, y_interp_[i], tol_[0]); + EXPECT_NEAR(dytmp, dy_interp_[i], tol_[1]); + EXPECT_NEAR(d2ytmp, d2y_interp_[i], tol_[2]); + } } } -TEST_F(CubicSplineTest, NotAKnot_sinx_Notuniform) -{ - CubicSpline cubspl; - int n = 10; - x_[0] = 0; - y_[0] = sin(x_[0]); - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i) * 2 * PI; - y_[9 - i] = sin(x_[9 - i]); - } - auto f = [](double x) -> double { return sin(x); }; - double err = count_err(x_, f, n); - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 3; - x_interp_[4] = 4.559; - x_interp_[5] = 6; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.009999833334167; - y_ref_[2] = 0.913413361341225; - y_ref_[3] = 0.141120008059867; - y_ref_[4] = -0.988258957900347; - y_ref_[5] = -0.279415498198926; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, FirstDeriv_2x_uniform) +TEST_F(CubicSplineTest, ErrorBound) { - CubicSpline cubspl; + // Error bound formula used in this test correspond to the complete cubic + // spline interpolant (exact first_deriv boundary conditions at both ends). + // This does not apply to other types of boundary conditions. - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = ((double)i * 10.0); - y_[i] = 2 * x_[i]; - } - auto f = [](double x) -> double { return 0; }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - 2, - 2); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 1.60; - x_interp_[2] = 3.20; - x_interp_[3] = 4.80; - x_interp_[4] = 6.40; - x_interp_[5] = 8.00; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 3.200000000000000; - y_ref_[2] = 6.400000000000000; - y_ref_[3] = 9.600000000000000; - y_ref_[4] = 12.800000000000000; - y_ref_[5] = 16.000000000000000; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} + // knots (logspace) + int n = 100; + double xmin = 0.1; + double xmax = 10; -TEST_F(CubicSplineTest, FirstDeriv_2x_Notuniform) -{ - CubicSpline cubspl; + double rho0 = std::log(xmin); + double drho = (std::log(xmax) - rho0) / (n - 1); + std::for_each(x_, x_ + n, [&](double& x) { x = std::exp(rho0 + (&x - x_) * drho); }); - int n = 10; - x_[0] = 0; - y_[0] = 0; - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i); - y_[9 - i] = 2 * (x_[9 - i]); - } - auto f = [](double x) -> double { return 0; }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - 2, - 2); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 0.10; - x_interp_[3] = 0.33; - x_interp_[4] = 0.66; - x_interp_[5] = 0.99; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.020000000000000; - y_ref_[2] = 0.200000000000000; - y_ref_[3] = 0.660000000000000; - y_ref_[4] = 1.320000000000000; - y_ref_[5] = 1.980000000000000; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} + // places to evaluate the interpolant + int n_interp = 777; + double dx_interp = (xmax - xmin) / (n_interp - 1); + std::for_each(x_interp_, x_interp_ + n_interp, + [&](double& x) { x = (&x - x_interp_) * dx_interp + xmin; }); -TEST_F(CubicSplineTest, SecondDeriv_2x_uniform) -{ - CubicSpline cubspl; + // make sure x_interp is inside the range of x + x_interp_[0] += tol_[0]; + x_interp_[n_interp - 1] -= tol_[0]; - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = ((double)i * 10.0); - y_[i] = 2 * x_[i]; - } - auto f = [](double x) -> double { return 0; }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - 0, - 0); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 1.60; - x_interp_[2] = 3.20; - x_interp_[3] = 4.80; - x_interp_[4] = 6.40; - x_interp_[5] = 8.00; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 3.200000000000000; - y_ref_[2] = 6.400000000000000; - y_ref_[3] = 9.600000000000000; - y_ref_[4] = 12.800000000000000; - y_ref_[5] = 16.000000000000000; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) + for (size_t i = 0; i < f_.size(); ++i) { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} + std::transform(x_, x_ + n, y_, f_[i][0]); -TEST_F(CubicSplineTest, SecondDeriv_2x_Notuniform) -{ - CubicSpline cubspl; + // complete cubic spline (exact first_deriv boundary conditions at both ends) + CubicSpline::build( + n, x_, y_, + {BoundaryType::first_deriv, f_[i][1](x_[0])}, + {BoundaryType::first_deriv, f_[i][1](x_[n - 1])}, + dy_ + ); - int n = 10; - x_[0] = 0; - y_[0] = 0; - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i); - y_[9 - i] = 2 * (x_[9 - i]); - } - auto f = [](double x) -> double { return 0; }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - 0, - 0); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 0.10; - x_interp_[3] = 0.33; - x_interp_[4] = 0.66; - x_interp_[5] = 0.99; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.020000000000000; - y_ref_[2] = 0.200000000000000; - y_ref_[3] = 0.660000000000000; - y_ref_[4] = 1.320000000000000; - y_ref_[5] = 1.980000000000000; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} + CubicSpline::eval( + n, x_, y_, dy_, + n_interp, x_interp_, y_interp_, dy_interp_, d2y_interp_ + ); -TEST_F(CubicSplineTest, NotAKnot_2x_uniform) -{ - CubicSpline cubspl; + double* diff[3] = {y_interp_, dy_interp_, d2y_interp_}; + for (int d = 0; d < 3; ++d) + { + std::transform(x_interp_, x_interp_ + n_interp, diff[d], diff[d], + [&](double x, double y) { return std::abs(y - f_[i][d](x)); }); - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = ((double)i * 10.0); - y_[i] = 2 * (x_[i]); - } - auto f = [](double x) -> double { return 0; }; - double err = count_err(x_, f, n); - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 1.60; - x_interp_[2] = 3.20; - x_interp_[3] = 4.80; - x_interp_[4] = 6.40; - x_interp_[5] = 8.00; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 3.200000000000000; - y_ref_[2] = 6.400000000000000; - y_ref_[3] = 9.600000000000000; - y_ref_[4] = 12.800000000000000; - y_ref_[5] = 16.000000000000000; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); + double err_bound = error_bound(n, x_, f_[i][4], d); + EXPECT_TRUE(std::all_of(diff[d], diff[d] + n_interp, + [err_bound](double diff) { return diff < err_bound; })); + } } } -TEST_F(CubicSplineTest, NotAKnot_2x_Notuniform) -{ - CubicSpline cubspl; - - int n = 10; - x_[0] = 0; - y_[0] = 0; - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i); - y_[9 - i] = 2 * (x_[9 - i]); - } - - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - auto f = [](double x) -> double { return 0; }; - double err = count_err(x_, f, n); - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 0.10; - x_interp_[3] = 0.33; - x_interp_[4] = 0.66; - x_interp_[5] = 0.99; - - y_ref_[0] = 0.000000000000000; - y_ref_[1] = 0.020000000000000; - y_ref_[2] = 0.200000000000000; - y_ref_[3] = 0.660000000000000; - y_ref_[4] = 1.320000000000000; - y_ref_[5] = 1.980000000000000; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} -TEST_F(CubicSplineTest, FirstDeriv_expx_uniform) +TEST_F(CubicSplineTest, Reserve) { - CubicSpline cubspl; - - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = ((double)i / 10.0); - y_[i] = exp(-x_[i]); - } - auto f = [](double x) -> double { return exp(-x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - -1, - -0.406570); - - int ni = 6; - x_interp_[0] = 0.000000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 0.200000000000000; - x_interp_[4] = 0.400000000000000; - x_interp_[5] = 0.800000000000000; - y_ref_[0] = 1.000000000000000; - y_ref_[1] = 0.990049833749168; - y_ref_[2] = 0.904837418035960; - y_ref_[3] = 0.818730753077982; - y_ref_[4] = 0.670320046035639; - y_ref_[5] = 0.449328964117222; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) + int n_spline = 20; + int n = 1000; + double x0 = 0.0, dx = 0.01; + for (int i = 0; i < n; ++i) { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); + x_[i] = x0 + i * dx; + y_[i] = std::sin(x_[i]); } -} -TEST_F(CubicSplineTest, FirstDeriv_expx_Notuniform) -{ - CubicSpline cubspl; + CubicSpline cubspl(n, x0, dx, y_); + cubspl.reserve(n_spline); + EXPECT_EQ(cubspl.heap_usage(), n_spline * 2 * n * sizeof(double)); - int n = 10; - x_[0] = 0; - y_[0] = 1; - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i); - y_[9 - i] = exp(-(x_[9 - i])); - } - auto f = [](double x) -> double { return exp(-x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - -1, - -0.367879); - - int ni = 6; - - x_interp_[0] = 0.000000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 0.330000000000000; - x_interp_[4] = 0.660000000000000; - x_interp_[5] = 0.990000000000000; - y_ref_[0] = 1.000000000000000; - y_ref_[1] = 0.990049833749168; - y_ref_[2] = 0.904837418035960; - y_ref_[3] = 0.718923733431926; - y_ref_[4] = 0.516851334491699; - y_ref_[5] = 0.371576691022046; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } + cubspl = CubicSpline(n, x_, y_); + cubspl.reserve(n_spline); + EXPECT_EQ(cubspl.heap_usage(), (1 + n_spline * 2) * n * sizeof(double)); } -TEST_F(CubicSplineTest, SecondDeriv_expx_uniform) -{ - CubicSpline cubspl; - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = ((double)i / 10.0); - y_[i] = exp(-x_[i]); - } - auto f = [](double x) -> double { return exp(-x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - 1, - 0.406570); - - int ni = 6; - x_interp_[0] = 0.000000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 0.200000000000000; - x_interp_[4] = 0.400000000000000; - x_interp_[5] = 0.800000000000000; - y_ref_[0] = 1.000000000000000; - y_ref_[1] = 0.990049833749168; - y_ref_[2] = 0.904837418035960; - y_ref_[3] = 0.818730753077982; - y_ref_[4] = 0.670320046035639; - y_ref_[5] = 0.449328964117222; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, SecondDeriv_expx_Notuniform) +TEST_F(CubicSplineTest, MinMax) { - CubicSpline cubspl; - int n = 10; - x_[0] = 0; - y_[0] = 1; - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i); - y_[9 - i] = exp(-(x_[9 - i])); - } - auto f = [](double x) -> double { return exp(-x); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - 1, - 0.367879441171442); - - int ni = 6; - x_interp_[0] = 0.000000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 0.200000000000000; - x_interp_[4] = 0.400000000000000; - x_interp_[5] = 0.800000000000000; - - y_ref_[0] = 1.000000000000000; - y_ref_[1] = 0.990049833749168; - y_ref_[2] = 0.904837418035960; - y_ref_[3] = 0.818730753077982; - y_ref_[4] = 0.670320046035639; - y_ref_[5] = 0.449328964117222; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) + int n = 1000; + double x0 = 0.0, dx = 0.01; + for (int i = 0; i < n; ++i) { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); + x_[i] = x0 + i * dx; + y_[i] = std::sin(x_[i]); } -} -// TEST_F(CubicSplineTest, NotAKnot_expx_uniform) -// { -// CubicSpline cubspl; - -// int n = 10; -// for (int i = 0; i != n; ++i) -// { -// x_[i] = ((double)i / 10.0); -// y_[i] = exp(-x_[i]); -// } -// double err = count_err(x_, 2, n); - -// cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - -// int ni = 6; - -// x_interp_[0] = 0.000000000000000; -// x_interp_[1] = 0.010000000000000; -// x_interp_[2] = 0.100000000000000; -// x_interp_[3] = 0.200000000000000; -// x_interp_[4] = 0.400000000000000; -// x_interp_[5] = 0.800000000000000; -// y_ref_[0] = 1.000000000000000; -// y_ref_[1] = 0.990049833749168; -// y_ref_[2] = 0.904837418035960; -// y_ref_[3] = 0.818730753077982; -// y_ref_[4] = 0.670320046035639; -// y_ref_[5] = 0.449328964117222; - -// cubspl.eval(ni, x_interp_, y_interp_); -// for (int i = 0; i != ni; ++i) -// { -// EXPECT_NEAR(y_interp_[i], y_ref_[i], err); -// } -// } - -TEST_F(CubicSplineTest, NotAKnot_expx_Notuniform) -{ - CubicSpline cubspl; + CubicSpline cubspl(n, x_, y_); + EXPECT_EQ(cubspl.xmin(), x_[0]); + EXPECT_EQ(cubspl.xmax(), x_[n - 1]); - int n = 10; - x_[0] = 0; - y_[0] = 1; - for (int i = 8; i >= 0; i--) - { - x_[9 - i] = exp(-i); - y_[9 - i] = exp(-(x_[9 - i])); - } - auto f = [](double x) -> double { return exp(-x); }; - double err = count_err(x_, f, n); - - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - - int ni = 6; - x_interp_[0] = 0.000000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 0.200000000000000; - x_interp_[4] = 0.400000000000000; - x_interp_[5] = 0.800000000000000; - y_ref_[0] = 1.000000000000000; - y_ref_[1] = 0.990049833749168; - y_ref_[2] = 0.904837418035960; - y_ref_[3] = 0.818730753077982; - y_ref_[4] = 0.670320046035639; - y_ref_[5] = 0.449328964117222; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, FirstDeriv_lnx_uniform) -{ - CubicSpline cubspl; + int m = 300; + cubspl = CubicSpline(m, x0, dx, y_); + EXPECT_EQ(cubspl.xmin(), x0); + EXPECT_EQ(cubspl.xmax(), x0 + (m - 1) * dx); - int n = 10; - for (int i = 1; i != n + 1; ++i) - { - x_[i] = ((double)i / 2.0); - y_[i] = log(x_[i]); - } - auto f = [](double x) -> double { return -6 * pow(x, -4); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - 2, - 0.2); - - int ni = 6; - x_interp_[0] = 0.001000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 1.000000000000000; - x_interp_[4] = 2.000000000000000; - x_interp_[5] = 4.000000000000000; - y_ref_[0] = -6.907755278982137; - y_ref_[1] = -4.605170185988091; - y_ref_[2] = -2.302585092994045; - y_ref_[3] = 0.000000000000000; - y_ref_[4] = 0.693147180559945; - y_ref_[5] = 1.386294361119891; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } } -TEST_F(CubicSplineTest, FirstDeriv_lnx_Notuniform) -{ - CubicSpline cubspl; - int n = 10; - for (int i = 9; i >= 0; i--) - { - x_[9 - i] = exp(-i) * 5; - y_[9 - i] = log(x_[9 - i]); - } - auto f = [](double x) -> double { return -6 * pow(x, -4); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::first_deriv, - CubicSpline::BoundaryCondition::first_deriv, - 1620.616785515076799, - 0.2); - - int ni = 6; - x_interp_[0] = 0.001000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 1.000000000000000; - x_interp_[4] = 2.000000000000000; - x_interp_[5] = 4.000000000000000; - y_ref_[0] = -6.907755278982137; - y_ref_[1] = -4.605170185988091; - y_ref_[2] = -2.302585092994045; - y_ref_[3] = 0.000000000000000; - y_ref_[4] = 0.693147180559945; - y_ref_[5] = 1.386294361119891; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, SecondDeriv_lnx_uniform) +TEST_F(CubicSplineTest, CrossCheck) { - CubicSpline cubspl; + std::vector fnames = { + "./data/sin_not_a_knot.dat", + "./data/cos_periodic.dat", + "./data/exp_first_deriv.dat", + "./data/log_second_deriv.dat", + "./data/sqrt_mix_bc.dat", + "./data/two_points_periodic.dat", + "./data/two_points_first_deriv.dat", + "./data/two_points_second_deriv.dat", + "./data/three_points_not_a_knot.dat", + }; - int n = 10; - for (int i = 1; i != n + 1; ++i) - { - x_[i] = ((double)i / 2.0); - y_[i] = log(x_[i]); - } - auto f = [](double x) -> double { return -6 * pow(x, -4); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - -4, - -0.04); - - int ni = 6; - x_interp_[0] = 0.001000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 1.000000000000000; - x_interp_[4] = 2.000000000000000; - x_interp_[5] = 4.000000000000000; - y_ref_[0] = -6.907755278982137; - y_ref_[1] = -4.605170185988091; - y_ref_[2] = -2.302585092994045; - y_ref_[3] = 0.000000000000000; - y_ref_[4] = 0.693147180559945; - y_ref_[5] = 1.386294361119891; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} + int n = 0, n_interp = 0; + BoundaryCondition bc_start, bc_end; -TEST_F(CubicSplineTest, SecondDeriv_lnx_Notuniform) -{ - CubicSpline cubspl; - int n = 10; - for (int i = 9; i >= 0; i--) - { - x_[9 - i] = exp(-i) * 5; - y_[9 - i] = log(x_[9 - i]); - } - auto f = [](double x) -> double { return -6 * pow(x, -4); }; - double err = count_err(x_, f, n); - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - -2626398.765493220183998, - 0.04); - - int ni = 6; - x_interp_[0] = 0.001000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 1.000000000000000; - x_interp_[4] = 2.000000000000000; - x_interp_[5] = 4.000000000000000; - y_ref_[0] = -6.907755278982137; - y_ref_[1] = -4.605170185988091; - y_ref_[2] = -2.302585092994045; - y_ref_[3] = 0.000000000000000; - y_ref_[4] = 0.693147180559945; - y_ref_[5] = 1.386294361119891; - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) + for (const auto& fname : fnames) { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, NotAKnot_lnx_uniform) -{ - CubicSpline cubspl; + read(fname, n, x_, y_, bc_start, bc_end, + n_interp, x_interp_, y_ref_, dy_ref_, d2y_ref_); + CubicSpline cubspl(n, x_, y_, bc_start, bc_end); + cubspl.eval(n_interp, x_interp_, y_interp_, dy_interp_, d2y_interp_); - int n = 10; - for (int i = 1; i != n + 1; ++i) - { - x_[i] = ((double)i / 2.0); - y_[i] = log(x_[i]); - } - auto f = [](double x) -> double { return -6 * pow(x, -4); }; - double err = count_err(x_, f, n); - - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - - int ni = 6; - x_interp_[0] = 0.001000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 1.000000000000000; - x_interp_[4] = 2.000000000000000; - x_interp_[5] = 4.000000000000000; - y_ref_[0] = -6.907755278982137; - y_ref_[1] = -4.605170185988091; - y_ref_[2] = -2.302585092994045; - y_ref_[3] = 0.000000000000000; - y_ref_[4] = 0.693147180559945; - y_ref_[5] = 1.386294361119891; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); + double* diff[] = {y_interp_, dy_interp_, d2y_interp_}; + double* ref[] = {y_ref_, dy_ref_, d2y_ref_}; + for (int d = 0; d < 3; ++d) + { + std::transform(diff[d], diff[d] + n_interp, ref[d], diff[d], std::minus()); + EXPECT_TRUE(std::all_of(diff[d], diff[d] + n_interp, + [this, d](double diff) { return std::abs(diff) < tol_[d]; })); + } } } -TEST_F(CubicSplineTest, NotAKnot_lnx_Notuniform) -{ - CubicSpline cubspl; - int n = 10; - for (int i = 9; i >= 0; i--) - { - x_[9 - i] = exp(-i) * 5; - y_[9 - i] = log(x_[9 - i]); - } - auto f = [](double x) -> double { return -6 * pow(x, -4); }; - double err = count_err(x_, f, n); - - cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); - - int ni = 6; - x_interp_[0] = 0.001000000000000; - x_interp_[1] = 0.010000000000000; - x_interp_[2] = 0.100000000000000; - x_interp_[3] = 1.000000000000000; - x_interp_[4] = 2.000000000000000; - x_interp_[5] = 4.000000000000000; - y_ref_[0] = -6.907755278982137; - y_ref_[1] = -4.605170185988091; - y_ref_[2] = -2.302585092994045; - y_ref_[3] = 0.000000000000000; - y_ref_[4] = 0.693147180559945; - y_ref_[5] = 1.386294361119891; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], err); - } -} - -TEST_F(CubicSplineTest, SecondDeriv) +int main() { - CubicSpline cubspl; - - int n = 10; - for (int i = 0; i != n; ++i) - { - x_[i] = std::sqrt(i); - y_[i] = std::exp(-x_[i]); - } - - cubspl.build(n, - x_, - y_, - CubicSpline::BoundaryCondition::second_deriv, - CubicSpline::BoundaryCondition::second_deriv, - -1, - -0.01); - - int ni = 6; - x_interp_[0] = 0.0; - x_interp_[1] = 0.01; - x_interp_[2] = 1.99; - x_interp_[3] = 2.0; - x_interp_[4] = 2.54; - x_interp_[5] = 3.00; - - y_ref_[0] = 1.; - y_ref_[1] = 0.9952111318752899; - y_ref_[2] = 0.13668283949289; - y_ref_[3] = 0.1353352832366127; - y_ref_[4] = 0.0788753653329337; - y_ref_[5] = 0.0497870683678639; - - cubspl.eval(ni, x_interp_, y_interp_); - for (int i = 0; i != ni; ++i) - { - EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); - } + ::testing::InitGoogleTest(); + return RUN_ALL_TESTS(); } -int main(int argc, char** argv) -{ - -#ifdef __MPI - MPI_Init(&argc, &argv); - MPI_Comm_size(MPI_COMM_WORLD, &GlobalV::NPROC); - MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK); -#endif - testing::InitGoogleTest(&argc, argv); - int result = RUN_ALL_TESTS(); - -#ifdef __MPI - MPI_Finalize(); -#endif - - return result; -} diff --git a/source/module_base/test/data/cos_periodic.dat b/source/module_base/test/data/cos_periodic.dat new file mode 100644 index 0000000000..2d0ae6b062 --- /dev/null +++ b/source/module_base/test/data/cos_periodic.dat @@ -0,0 +1,7 @@ +periodic periodic + 0.00000000000000000e+00 1.00000000000000006e-01 1.00000000000000000e+00 1.50000000000000000e+00 3.00000000000000000e+00 6.28318530717958623e+00 + 1.00000000000000000e+00 9.95004165278025821e-01 5.40302305868139765e-01 7.07372016677029064e-02 -9.89992496600445415e-01 1.00000000000000000e+00 + 0.00000000000000000e+00 1.74532925199432948e-01 3.49065850398865896e-01 5.23598775598298816e-01 6.98131700797731791e-01 8.72664625997164767e-01 1.04719755119659763e+00 1.22173047639603061e+00 1.39626340159546358e+00 1.57079632679489656e+00 1.74532925199432953e+00 1.91986217719376251e+00 2.09439510239319526e+00 2.26892802759262846e+00 2.44346095279206121e+00 2.61799387799149441e+00 2.79252680319092716e+00 2.96705972839035992e+00 3.14159265358979312e+00 3.31612557878922587e+00 3.49065850398865907e+00 3.66519142918809182e+00 3.83972435438752502e+00 4.01425727958695777e+00 4.18879020478639053e+00 4.36332312998582328e+00 4.53785605518525692e+00 4.71238898038468967e+00 4.88692190558412243e+00 5.06145483078355518e+00 5.23598775598298882e+00 5.41052068118242158e+00 5.58505360638185433e+00 5.75958653158128708e+00 5.93411945678071984e+00 6.10865238198015348e+00 6.28318530717958623e+00 + 1.00000000000000000e+00 9.84381299886306516e-01 9.38016522020507670e-01 8.63713424115361450e-01 7.64165387715657451e-01 6.42065794366185383e-01 5.00118352188328785e-01 3.42056178626561569e-01 1.73503016854116870e-01 2.70592234453568648e-04 -1.72386293431320586e-01 -3.39925662651143046e-01 -4.97852121855354546e-01 -6.41670277474296391e-01 -7.66884735938309081e-01 -8.69000103677733948e-01 -9.43520987122911547e-01 -9.85951992704182878e-01 -9.92542393974230608e-01 -9.64850207988183883e-01 -9.06738697819773454e-01 -8.22080502627916765e-01 -7.14748261571530152e-01 -5.88614613809530951e-01 -4.47552198500836496e-01 -2.95433654804362900e-01 -1.36131621879026887e-01 2.64812611162533740e-02 1.88532355022561493e-01 3.46149020680980524e-01 4.95458618932594463e-01 6.32588510618484978e-01 7.53666056579735066e-01 8.54818617657429058e-01 9.32173554692648398e-01 9.81858228526478083e-01 1.00000000000000000e+00 + 6.18473891098361445e-03 -1.80469237892337431e-01 -3.48259990143118536e-01 -5.00618802880412384e-01 -6.37545676104219083e-01 -7.59040609814538803e-01 -8.64447219749984619e-01 -9.41247940672875094e-01 -9.84665089228650259e-01 -9.95326120656199964e-01 -9.78883265201603403e-01 -9.36683701867835450e-01 -8.68727430654895882e-01 -7.75014451562784812e-01 -6.55544764591502682e-01 -5.10318369741048716e-01 -3.39335267011423580e-01 -1.42595456402626830e-01 6.41233977052624166e-02 2.49507580803381257e-01 4.12703175963197166e-01 5.53710183184709215e-01 6.72528602467918124e-01 7.69158433812823450e-01 8.43599677219425303e-01 8.95852332687724018e-01 9.25916400217719260e-01 9.33791879809410919e-01 9.19478771462799216e-01 8.82977075177883819e-01 8.24286790954665616e-01 7.43407918793143718e-01 6.40340458693318570e-01 5.15084410655190172e-01 3.67639774678758080e-01 1.98006550764022293e-01 6.18473891098361445e-03 +-1.16262364588415146e+00 -1.00557944472069671e+00 -9.17160944338298312e-01 -8.28742443955899910e-01 -7.40323943573501397e-01 -6.51905443191102996e-01 -5.35672607329632511e-01 -3.44398827159106469e-01 -1.53125046988580316e-01 2.04230892878108500e-02 1.67998154850598658e-01 3.15573220413386479e-01 4.63148285976174190e-01 6.10723351538962067e-01 7.58298417101749722e-01 9.05873482664537821e-01 1.05344854822732525e+00 1.20102361379011291e+00 1.12573875011125457e+00 9.98607505889287705e-01 8.71476261667320506e-01 7.44345017445353752e-01 6.17213773223386442e-01 4.90082529001419687e-01 3.62951284779452932e-01 2.35820040557485955e-01 1.08688796335518534e-01 -1.84424478864484431e-02 -1.45573692108415198e-01 -2.72704936330382175e-01 -3.99836180552349374e-01 -5.26967424774316351e-01 -6.54098668996283328e-01 -7.81229913218249861e-01 -9.08361157440217282e-01 -1.03549240166218426e+00 -1.16262364588415146e+00 diff --git a/source/module_base/test/data/exp_first_deriv.dat b/source/module_base/test/data/exp_first_deriv.dat new file mode 100644 index 0000000000..e5144b72b0 --- /dev/null +++ b/source/module_base/test/data/exp_first_deriv.dat @@ -0,0 +1,7 @@ +first_deriv:1.10517091807564771e+00 first_deriv:2.71828182845904509e+00 + 1.00000000000000006e-01 1.29154966501488389e-01 1.66810053720005874e-01 2.15443469003188337e-01 2.78255940220712428e-01 3.59381366380462752e-01 4.64158883361277863e-01 5.99484250318940926e-01 7.74263682681126997e-01 1.00000000000000000e+00 + 1.10517091807564771e+00 1.13786644168535322e+00 1.18152981679905356e+00 1.24041185922219421e+00 1.32082420599103290e+00 1.43244298301307027e+00 1.59067568226275657e+00 1.82117928549843455e+00 2.16899447087117636e+00 2.71828182845904509e+00 + 1.00000000000000006e-01 1.25000000000000000e-01 1.50000000000000022e-01 1.75000000000000017e-01 2.00000000000000011e-01 2.25000000000000006e-01 2.50000000000000000e-01 2.75000000000000022e-01 3.00000000000000044e-01 3.25000000000000011e-01 3.49999999999999978e-01 3.75000000000000000e-01 4.00000000000000022e-01 4.25000000000000044e-01 4.50000000000000067e-01 4.74999999999999978e-01 5.00000000000000000e-01 5.25000000000000022e-01 5.50000000000000044e-01 5.75000000000000067e-01 5.99999999999999978e-01 6.25000000000000000e-01 6.50000000000000022e-01 6.75000000000000044e-01 7.00000000000000067e-01 7.24999999999999978e-01 7.50000000000000000e-01 7.75000000000000022e-01 8.00000000000000044e-01 8.25000000000000067e-01 8.49999999999999978e-01 8.75000000000000000e-01 9.00000000000000022e-01 9.25000000000000044e-01 9.50000000000000067e-01 9.74999999999999978e-01 1.00000000000000000e+00 + 1.10517091807564771e+00 1.13314845305620615e+00 1.16183423781052575e+00 1.19124621006218145e+00 1.22140274895196677e+00 1.25232269983173161e+00 1.28402537693915098e+00 1.31653067826587877e+00 1.34985870568413047e+00 1.38403052288246253e+00 1.41906754700129101e+00 1.45499127838897557e+00 1.49182434073154035e+00 1.52959016783122848e+00 1.56831221140395694e+00 1.60801396407294228e+00 1.64872023299870918e+00 1.69045738865219253e+00 1.73325189263525137e+00 1.77713020654974518e+00 1.82211879200369298e+00 1.86824485646834759e+00 1.91553841190619289e+00 1.96403012931547272e+00 2.01375067969443089e+00 2.06473073404130947e+00 2.11700096335435450e+00 2.17059203866204742e+00 2.22553592214226725e+00 2.28186930422749290e+00 2.33962995740218638e+00 2.39885565415081015e+00 2.45958416695782534e+00 2.52185326830769485e+00 2.58570073068487982e+00 2.65116432657384271e+00 2.71828182845904465e+00 + 1.10517091807564771e+00 1.13314857626511545e+00 1.16183448755634799e+00 1.19124518701385762e+00 1.22140376069562406e+00 1.25232046400808739e+00 1.28402663663421768e+00 1.31653040155273815e+00 1.34985411865053440e+00 1.38403404177963485e+00 1.41907067229022621e+00 1.45497988810568324e+00 1.49182144407227635e+00 1.52960106867554990e+00 1.56831876191550346e+00 1.60798584383446830e+00 1.64869208090593933e+00 1.69045678205178129e+00 1.73327994727199441e+00 1.77716157656657869e+00 1.82210170576561081e+00 1.86818792472837214e+00 1.91550098362087740e+00 1.96404088244312680e+00 2.01380762119512058e+00 2.06480119987685740e+00 2.11702161848833859e+00 2.17046900023482614e+00 2.22529349474582894e+00 2.28162888863530977e+00 2.33947518190326731e+00 2.39883237454970333e+00 2.45970046657461650e+00 2.52207945797800726e+00 2.58596934875987605e+00 2.65137013892022200e+00 2.71828182845904553e+00 + 1.10510282004738825e+00 1.13310983511001884e+00 1.16189187653084680e+00 1.19124256407428741e+00 1.22144333046701914e+00 1.25229505919741335e+00 1.28419875089301683e+00 1.31610244258862052e+00 1.35006277753420334e+00 1.38433107279383738e+00 1.41859936805347142e+00 1.45490086593012635e+00 1.49242361139733326e+00 1.52994635686453995e+00 1.56746910233174686e+00 1.60708020137141450e+00 1.64941876434625834e+00 1.69175732732110240e+00 1.73409589029594624e+00 1.77643445327079030e+00 1.81891195991557297e+00 1.86798555710533210e+00 1.91705915429509122e+00 1.96613275148485056e+00 2.01520634867460968e+00 2.06427994586436858e+00 2.11335354305412793e+00 2.16276179287056580e+00 2.22319776800966951e+00 2.28363374314877277e+00 2.34406971828787603e+00 2.40450569342697928e+00 2.46494166856608299e+00 2.52537764370518625e+00 2.58581361884428951e+00 2.64624959398339277e+00 2.70668556912249603e+00 diff --git a/source/module_base/test/data/gen_ref.py b/source/module_base/test/data/gen_ref.py new file mode 100644 index 0000000000..abc1969657 --- /dev/null +++ b/source/module_base/test/data/gen_ref.py @@ -0,0 +1,64 @@ +''' +This script generates reference data for the cross-check in cubic_spline_test.cpp. +To add more tests, append a new Case object to the cases list. + +''' +import numpy as np +from scipy.interpolate import CubicSpline + +class Case: + def __init__(self, f, x, bc, fname): + self.x = x + self.y = f(x) + self.bc = bc + self.fname = fname + if bc == 'periodic': + self.y[-1] = self.y[0] + +cases = [ + Case(np.sin, np.logspace(-1, 0, 10), 'not-a-knot', + 'sin_not_a_knot.dat'), + Case(np.cos, [0, 0.1, 1, 1.5, 3, 2*np.pi], 'periodic', + 'cos_periodic.dat'), + Case(np.exp, np.logspace(-1, 0, 10), ((1, np.exp(0.1)), (1, np.exp(1))), + 'exp_first_deriv.dat'), + Case(np.log, np.logspace(0, 1, 10), ((2, -1), (2, -0.01)), + 'log_second_deriv.dat'), + Case(np.sqrt, np.logspace(0, 1, 10), ((1, 0.5), 'not-a-knot'), + 'sqrt_mix_bc.dat'), + Case(lambda x: np.ones_like(x), [0, 1], 'periodic', + 'two_points_periodic.dat'), + Case(np.arccos, [0, 0.5], ((1, -1), (1, -2/np.sqrt(3))), + 'two_points_first_deriv.dat'), + Case(np.sqrt, [1, 4], ((2, 0.5), (2, 0.25)), + 'two_points_second_deriv.dat'), + Case(np.sin, [0, 0.3, 1], 'not-a-knot', + 'three_points_not_a_knot.dat'), +] + +n_interp = 37 + +for case in cases: + cubspl = CubicSpline(case.x, case.y, bc_type=case.bc) + x_interp = np.linspace(case.x[0], case.x[-1], n_interp) + y_interp = cubspl(x_interp) + dy_interp = cubspl(x_interp, 1) + d2y_interp = cubspl(x_interp, 2) + + with open(case.fname, 'w') as f: + if case.bc == 'periodic' or case.bc == 'not-a-knot': + f.write('{bc} {bc}\n'.format(bc=case.bc)) + else: + for b in case.bc: + if b == 'not-a-knot': + f.write('not-a-knot ') + else: + tag = 'first_deriv' if b[0] == 1 else 'second_deriv' + f.write('{tag}:{val:<24.17e} '.format(tag=tag, val=b[1])) + f.write('\n') + + for data in [case.x, case.y, x_interp, y_interp, dy_interp, d2y_interp]: + for elem in data: + f.write('% 24.17e '%(elem)) + f.write('\n') + diff --git a/source/module_base/test/data/log_second_deriv.dat b/source/module_base/test/data/log_second_deriv.dat new file mode 100644 index 0000000000..62060d5f77 --- /dev/null +++ b/source/module_base/test/data/log_second_deriv.dat @@ -0,0 +1,7 @@ +second_deriv:-1.00000000000000000e+00 second_deriv:-1.00000000000000002e-02 + 1.00000000000000000e+00 1.29154966501488389e+00 1.66810053720005880e+00 2.15443469003188381e+00 2.78255940220712450e+00 3.59381366380462763e+00 4.64158883361277841e+00 5.99484250318940859e+00 7.74263682681126930e+00 1.00000000000000000e+01 + 0.00000000000000000e+00 2.55842788110449526e-01 5.11685576220899052e-01 7.67528364331348634e-01 1.02337115244179810e+00 1.27921394055224780e+00 1.53505672866269705e+00 1.79089951677314629e+00 2.04674230488359621e+00 2.30258509299404590e+00 + 1.00000000000000000e+00 1.25000000000000000e+00 1.50000000000000000e+00 1.75000000000000000e+00 2.00000000000000000e+00 2.25000000000000000e+00 2.50000000000000000e+00 2.75000000000000000e+00 3.00000000000000000e+00 3.25000000000000000e+00 3.50000000000000000e+00 3.75000000000000000e+00 4.00000000000000000e+00 4.25000000000000000e+00 4.50000000000000000e+00 4.75000000000000000e+00 5.00000000000000000e+00 5.25000000000000000e+00 5.50000000000000000e+00 5.75000000000000000e+00 6.00000000000000000e+00 6.25000000000000000e+00 6.50000000000000000e+00 6.75000000000000000e+00 7.00000000000000000e+00 7.25000000000000000e+00 7.50000000000000000e+00 7.75000000000000000e+00 8.00000000000000000e+00 8.25000000000000000e+00 8.50000000000000000e+00 8.75000000000000000e+00 9.00000000000000000e+00 9.25000000000000000e+00 9.50000000000000000e+00 9.75000000000000000e+00 1.00000000000000000e+01 + 0.00000000000000000e+00 2.23201685990528503e-01 4.05482161721055334e-01 5.59647381157334922e-01 6.93212288502453533e-01 8.10938183318088091e-01 9.16348065656249267e-01 1.01160453727708877e+00 1.09865227971226798e+00 1.17872100366669685e+00 1.25277767969836296e+00 1.32176777395878364e+00 1.38635002079692837e+00 1.44697637249515854e+00 1.50409420887274625e+00 1.55814764821316754e+00 1.60947600084362952e+00 1.65829393447507156e+00 1.70480885095700785e+00 1.74922815213895189e+00 1.79175923970966466e+00 1.83259005049654222e+00 1.87183533226151133e+00 1.90959263384970646e+00 1.94595950410626362e+00 1.98103349187631750e+00 2.01491214600500346e+00 2.04769301510038781e+00 2.07946352550178926e+00 2.11027403528676105e+00 2.14016641953235709e+00 2.16918255331563126e+00 2.19736431171363611e+00 2.22475356980342642e+00 2.25139220266205564e+00 2.27732208536657810e+00 2.30258509299404635e+00 + 1.00254261613056794e+00 7.98334999625206598e-01 6.66752209256807338e-01 5.71952135021786212e-01 4.99597454958138398e-01 4.44699888470842797e-01 3.99955939076224509e-01 3.63472602732268946e-01 3.33587746439520605e-01 3.07606422584049077e-01 2.85491363057420178e-01 2.66837284058103641e-01 2.50118943859901521e-01 2.35190122938787882e-01 2.22050821294762640e-01 2.10610784319237826e-01 2.00154304624132307e-01 1.90527432327081331e-01 1.81730167428084816e-01 1.73762509927142816e-01 1.66624366318174705e-01 1.60087152540269423e-01 1.53940134142904944e-01 1.48183311126081213e-01 1.42816683489798313e-01 1.37840251234056188e-01 1.33254014358854894e-01 1.29057876274537020e-01 1.25134123654711077e-01 1.21377871343099950e-01 1.17789119339703638e-01 1.14367867644522142e-01 1.11114116257555462e-01 1.08027865178803598e-01 1.05109114408266549e-01 1.02357863945944316e-01 9.97741137918368992e-02 +-1.00000000000000377e+00 -6.33660932042887648e-01 -4.49235132334255161e-01 -3.25782694882304436e-01 -2.53054745626878019e-01 -1.95497023679798904e-01 -1.62454571477147702e-01 -1.29412119274496529e-01 -1.11657824079571477e-01 -9.61927667642009099e-02 -8.07277094488303704e-02 -7.04523993469851478e-02 -6.32943222386314758e-02 -5.61362451302777968e-02 -4.89781680219241247e-02 -4.34851335765309294e-02 -4.01667039843129556e-02 -3.68482743920949818e-02 -3.35298447998770149e-02 -3.02114152076590411e-02 -2.69292458727027513e-02 -2.53684643505395734e-02 -2.38076828283763989e-02 -2.22469013062132209e-02 -2.06861197840500430e-02 -1.91253382618868650e-02 -1.75645567397236871e-02 -1.60300110957334006e-02 -1.53600098628741376e-02 -1.46900086300148763e-02 -1.40200073971556132e-02 -1.33500061642963519e-02 -1.26800049314370888e-02 -1.20100036985778258e-02 -1.13400024657185645e-02 -1.06700012328593032e-02 -1.00000000000000401e-02 diff --git a/source/module_base/test/data/sin_not_a_knot.dat b/source/module_base/test/data/sin_not_a_knot.dat new file mode 100644 index 0000000000..417d0140ae --- /dev/null +++ b/source/module_base/test/data/sin_not_a_knot.dat @@ -0,0 +1,7 @@ +not-a-knot not-a-knot + 1.00000000000000006e-01 1.29154966501488389e-01 1.66810053720005874e-01 2.15443469003188337e-01 2.78255940220712428e-01 3.59381366380462752e-01 4.64158883361277863e-01 5.99484250318940926e-01 7.74263682681126997e-01 1.00000000000000000e+00 + 9.98334166468281548e-02 1.28796193418703991e-01 1.66037531159675150e-01 2.13780666055298912e-01 2.74679090976688467e-01 3.51695188663471381e-01 4.47670834718957189e-01 5.64216731736942645e-01 6.99189845133651677e-01 8.41470984807896505e-01 + 1.00000000000000006e-01 1.25000000000000000e-01 1.50000000000000022e-01 1.75000000000000017e-01 2.00000000000000011e-01 2.25000000000000006e-01 2.50000000000000000e-01 2.75000000000000022e-01 3.00000000000000044e-01 3.25000000000000011e-01 3.49999999999999978e-01 3.75000000000000000e-01 4.00000000000000022e-01 4.25000000000000044e-01 4.50000000000000067e-01 4.74999999999999978e-01 5.00000000000000000e-01 5.25000000000000022e-01 5.50000000000000044e-01 5.75000000000000067e-01 5.99999999999999978e-01 6.25000000000000000e-01 6.50000000000000022e-01 6.75000000000000044e-01 7.00000000000000067e-01 7.24999999999999978e-01 7.50000000000000000e-01 7.75000000000000022e-01 8.00000000000000044e-01 8.25000000000000067e-01 8.49999999999999978e-01 8.75000000000000000e-01 9.00000000000000022e-01 9.25000000000000044e-01 9.50000000000000067e-01 9.74999999999999978e-01 1.00000000000000000e+00 + 9.98334166468281548e-02 1.24674734514712698e-01 1.49438130358688637e-01 1.74108137455688206e-01 1.98669331470162691e-01 2.23106356655036864e-01 2.47403944178155422e-01 2.71546935540866008e-01 2.95520208213219726e-01 3.19308806225840347e-01 3.42897834817652569e-01 3.66272425041630545e-01 3.89418057463319811e-01 4.12320464705872358e-01 4.34965384966057578e-01 4.57338563261367681e-01 4.79425963790102938e-01 5.01213811410824883e-01 5.22688346176908625e-01 5.43835808141729826e-01 5.64642437361874472e-01 5.85094862720930697e-01 6.05181175115716008e-01 6.24889809005897301e-01 6.44209198851141918e-01 6.63127779111117310e-01 6.81633984245490487e-01 6.99716248713929012e-01 7.17363006976099671e-01 7.34562693491669916e-01 7.51303742720306977e-01 7.67574589121678086e-01 7.83363667155450361e-01 7.98659411281291143e-01 8.13450255958867552e-01 8.27724635647847040e-01 8.41470984807896505e-01 + 9.95004923184505952e-01 9.92197390241733368e-01 9.88771161271814858e-01 9.84726484485526354e-01 9.80066691602079154e-01 9.74793416196813478e-01 9.68912581684626351e-01 9.62425723364175201e-01 9.55336713363020729e-01 9.47651829867646733e-01 9.39371159698343727e-01 9.30499643073893545e-01 9.21055871973040285e-01 9.11041628742959020e-01 9.00456913383649749e-01 8.89303613354681888e-01 8.77596695966640139e-01 8.65339380713611495e-01 8.52531667595595843e-01 8.39173556612593186e-01 8.25265066443230833e-01 8.10851858679044746e-01 7.95976029301547761e-01 7.80637578310739877e-01 7.64836505706621206e-01 7.48572811489191747e-01 7.31846495658451279e-01 7.14657558214399802e-01 6.97005999157037537e-01 6.78891818486364262e-01 6.60315016202380312e-01 6.41275592305085462e-01 6.21773546794479715e-01 6.01808879670563179e-01 5.81381590933335857e-01 5.60491680582797747e-01 5.39139148618948627e-01 +-9.99273971679864453e-02 -1.24675238253821624e-01 -1.49423079339714465e-01 -1.74110306895194078e-01 -1.98673123780586686e-01 -2.23112904322202821e-01 -2.47353856652764159e-01 -2.71594808983325553e-01 -2.95479606336378597e-01 -3.19311073293541137e-01 -3.43142540250703731e-01 -3.66341401449573223e-01 -3.89160286618692008e-01 -4.11979171787810738e-01 -4.34798056956929524e-01 -4.57268738221928150e-01 -4.79284652821407420e-01 -5.01300567420886689e-01 -5.23316482020365958e-01 -5.45332396619845228e-01 -5.67275878301223169e-01 -5.85780742833658996e-01 -6.04285607366094824e-01 -6.22790471898530651e-01 -6.41295336430966478e-01 -6.59800200963402195e-01 -6.78305065495838022e-01 -6.96809930028275515e-01 -7.15314794560709233e-01 -7.33819659093142951e-01 -7.52324523625576669e-01 -7.70829388158010387e-01 -7.89334252690444216e-01 -8.07839117222877934e-01 -8.26343981755311652e-01 -8.44848846287745370e-01 -8.63353710820179088e-01 diff --git a/source/module_base/test/data/sqrt_mix_bc.dat b/source/module_base/test/data/sqrt_mix_bc.dat new file mode 100644 index 0000000000..1a63a03293 --- /dev/null +++ b/source/module_base/test/data/sqrt_mix_bc.dat @@ -0,0 +1,7 @@ +first_deriv:5.00000000000000000e-01 not-a-knot + 1.00000000000000000e+00 1.29154966501488389e+00 1.66810053720005880e+00 2.15443469003188381e+00 2.78255940220712450e+00 3.59381366380462763e+00 4.64158883361277841e+00 5.99484250318940859e+00 7.74263682681126930e+00 1.00000000000000000e+01 + 1.00000000000000000e+00 1.13646366638572482e+00 1.29154966501488389e+00 1.46779926762206947e+00 1.66810053720005880e+00 1.89573565240637598e+00 2.15443469003188381e+00 2.44843674682222678e+00 2.78255940220712450e+00 3.16227766016837952e+00 + 1.00000000000000000e+00 1.25000000000000000e+00 1.50000000000000000e+00 1.75000000000000000e+00 2.00000000000000000e+00 2.25000000000000000e+00 2.50000000000000000e+00 2.75000000000000000e+00 3.00000000000000000e+00 3.25000000000000000e+00 3.50000000000000000e+00 3.75000000000000000e+00 4.00000000000000000e+00 4.25000000000000000e+00 4.50000000000000000e+00 4.75000000000000000e+00 5.00000000000000000e+00 5.25000000000000000e+00 5.50000000000000000e+00 5.75000000000000000e+00 6.00000000000000000e+00 6.25000000000000000e+00 6.50000000000000000e+00 6.75000000000000000e+00 7.00000000000000000e+00 7.25000000000000000e+00 7.50000000000000000e+00 7.75000000000000000e+00 8.00000000000000000e+00 8.25000000000000000e+00 8.50000000000000000e+00 8.75000000000000000e+00 9.00000000000000000e+00 9.25000000000000000e+00 9.50000000000000000e+00 9.75000000000000000e+00 1.00000000000000000e+01 + 1.00000000000000000e+00 1.11803728698643989e+00 1.22475682769484950e+00 1.32287978023881414e+00 1.41422464207946574e+00 1.50000375623387994e+00 1.58115482205509994e+00 1.65831336044919975e+00 1.73206149587978508e+00 1.80279292004163660e+00 1.87083207865261780e+00 1.93649703530829753e+00 2.00002057410239775e+00 2.06157397820173083e+00 2.12132717083545463e+00 2.17944902427163800e+00 2.23607463862589828e+00 2.29129895056073130e+00 2.34521455546809943e+00 2.39791404873996372e+00 2.44949002569315066e+00 2.50002598372984819e+00 2.54957121153853006e+00 2.59816695899999894e+00 2.64585447599505841e+00 2.69267501240451246e+00 2.73866981810916421e+00 2.78388014298981767e+00 2.82834723692727730e+00 2.87211234980234575e+00 2.91521673149582661e+00 2.95770163188852342e+00 2.99960830086124020e+00 3.04097798829478050e+00 3.08185194406994789e+00 3.12227141806754682e+00 3.16227766016837952e+00 + 5.00000000000000000e-01 4.47098909750842544e-01 4.08224810882382050e-01 3.78053779833561365e-01 3.53481662512829198e-01 3.33407818352801533e-01 3.16209958323798612e-01 3.01467598935841907e-01 2.88742679905178157e-01 2.77324939467649523e-01 2.67204555498214813e-01 2.58260988871481856e-01 2.50040603634092751e-01 2.42499911313340782e-01 2.35638911909225951e-01 2.29428522778511251e-01 2.23638122316879184e-01 2.18218103423094306e-01 2.13168466097156617e-01 2.08489210339066144e-01 2.04180292443896155e-01 2.00134871770223183e-01 1.96254450619764736e-01 1.92539028992520816e-01 1.88988606888491395e-01 1.85603184307676500e-01 1.82382761250076131e-01 1.79327337715690260e-01 1.76436913704518944e-01 1.73711489216562182e-01 1.71151064251819918e-01 1.68755638810292180e-01 1.66525212891979024e-01 1.64459786496880367e-01 1.62559359624996236e-01 1.60823932276326687e-01 1.59253504450871608e-01 +-2.45211727308500382e-01 -1.77996994684759102e-01 -1.37423541149741701e-01 -1.07607040737062801e-01 -8.89698978287946640e-02 -7.37024413981044041e-02 -6.38804388339192647e-02 -5.40584362697341184e-02 -4.82656746863025995e-02 -4.30762488139266458e-02 -3.78868229415506921e-02 -3.42409267828305136e-02 -3.15221551162820982e-02 -2.88033834497336862e-02 -2.60846117831852707e-02 -2.39023649822226771e-02 -2.24208387108338972e-02 -2.09393124394451208e-02 -1.94577861680563409e-02 -1.79762598966675645e-02 -1.65116817411209045e-02 -1.58516836482628341e-02 -1.51916855554047638e-02 -1.45316874625466934e-02 -1.38716893696886247e-02 -1.32116912768305544e-02 -1.25516931839724857e-02 -1.18916950911143095e-02 -1.12316969982561819e-02 -1.05716989053980526e-02 -9.91170081253992322e-03 -9.25170271968179560e-03 -8.59170462682366798e-03 -7.93170653396553862e-03 -7.27170844110741014e-03 -6.61171034824928165e-03 -5.95171225539115403e-03 diff --git a/source/module_base/test/data/three_points_not_a_knot.dat b/source/module_base/test/data/three_points_not_a_knot.dat new file mode 100644 index 0000000000..0aaf9351c1 --- /dev/null +++ b/source/module_base/test/data/three_points_not_a_knot.dat @@ -0,0 +1,7 @@ +not-a-knot not-a-knot + 0.00000000000000000e+00 2.99999999999999989e-01 1.00000000000000000e+00 + 0.00000000000000000e+00 2.95520206661339546e-01 8.41470984807896505e-01 + 0.00000000000000000e+00 2.77777777777777762e-02 5.55555555555555525e-02 8.33333333333333287e-02 1.11111111111111105e-01 1.38888888888888895e-01 1.66666666666666657e-01 1.94444444444444420e-01 2.22222222222222210e-01 2.50000000000000000e-01 2.77777777777777790e-01 3.05555555555555525e-01 3.33333333333333315e-01 3.61111111111111105e-01 3.88888888888888840e-01 4.16666666666666630e-01 4.44444444444444420e-01 4.72222222222222210e-01 5.00000000000000000e-01 5.27777777777777790e-01 5.55555555555555580e-01 5.83333333333333259e-01 6.11111111111111049e-01 6.38888888888888840e-01 6.66666666666666630e-01 6.94444444444444420e-01 7.22222222222222210e-01 7.50000000000000000e-01 7.77777777777777679e-01 8.05555555555555469e-01 8.33333333333333259e-01 8.61111111111111049e-01 8.88888888888888840e-01 9.16666666666666630e-01 9.44444444444444420e-01 9.72222222222222210e-01 1.00000000000000000e+00 + 0.00000000000000000e+00 2.89141774610951995e-02 5.75117844399816239e-02 8.57928209366592559e-02 1.13757286951128120e-01 1.41405182483388209e-01 1.68736507533439495e-01 1.95751262101282020e-01 2.22449446186915784e-01 2.48831059790340758e-01 2.74896102911556972e-01 3.00644575550564286e-01 3.26076477707362977e-01 3.51191809381952769e-01 3.75990570574333827e-01 4.00472761284506151e-01 4.24638381512469687e-01 4.48487431258224378e-01 4.72019910521770392e-01 4.95235819303107561e-01 5.18135157602235941e-01 5.40717925419155421e-01 5.62984122753866445e-01 5.84933749606368458e-01 6.06566805976661683e-01 6.27883291864746340e-01 6.48883207270621987e-01 6.69566552194288955e-01 6.89933326635747024e-01 7.09983530594996526e-01 7.29717164072037017e-01 7.49134227066868941e-01 7.68234719579492076e-01 7.87018641609906422e-01 8.05485993158111757e-01 8.23636774224108636e-01 8.41470984807896394e-01 + 1.04660865727918528e+00 1.03521211991966933e+00 1.02381558256015315e+00 1.01241904520063719e+00 1.00102250784112101e+00 9.89625970481604944e-01 9.78229433122088876e-01 9.66832895762572919e-01 9.55436358403056851e-01 9.44039821043540783e-01 9.32643283684024715e-01 9.21246746324508647e-01 9.09850208964992579e-01 8.98453671605476623e-01 8.87057134245960555e-01 8.75660596886444487e-01 8.64264059526928530e-01 8.52867522167412462e-01 8.41470984807896394e-01 8.30074447448380437e-01 8.18677910088864369e-01 8.07281372729348412e-01 7.95884835369832344e-01 7.84488298010316276e-01 7.73091760650800319e-01 7.61695223291284251e-01 7.50298685931768183e-01 7.38902148572252226e-01 7.27505611212736158e-01 7.16109073853220091e-01 7.04712536493704134e-01 6.93315999134188066e-01 6.81919461774672109e-01 6.70522924415156041e-01 6.59126387055640084e-01 6.47729849696123905e-01 6.36333312336607948e-01 +-4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942578168e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 -4.10275344942577003e-01 diff --git a/source/module_base/test/data/two_points_first_deriv.dat b/source/module_base/test/data/two_points_first_deriv.dat new file mode 100644 index 0000000000..7ec6981d9c --- /dev/null +++ b/source/module_base/test/data/two_points_first_deriv.dat @@ -0,0 +1,7 @@ +first_deriv:-1.00000000000000000e+00 first_deriv:-1.15470053837925168e+00 + 0.00000000000000000e+00 5.00000000000000000e-01 + 1.57079632679489656e+00 1.04719755119659785e+00 + 0.00000000000000000e+00 1.38888888888888881e-02 2.77777777777777762e-02 4.16666666666666644e-02 5.55555555555555525e-02 6.94444444444444475e-02 8.33333333333333287e-02 9.72222222222222099e-02 1.11111111111111105e-01 1.25000000000000000e-01 1.38888888888888895e-01 1.52777777777777762e-01 1.66666666666666657e-01 1.80555555555555552e-01 1.94444444444444420e-01 2.08333333333333315e-01 2.22222222222222210e-01 2.36111111111111105e-01 2.50000000000000000e-01 2.63888888888888895e-01 2.77777777777777790e-01 2.91666666666666630e-01 3.05555555555555525e-01 3.19444444444444420e-01 3.33333333333333315e-01 3.47222222222222210e-01 3.61111111111111105e-01 3.75000000000000000e-01 3.88888888888888840e-01 4.02777777777777735e-01 4.16666666666666630e-01 4.30555555555555525e-01 4.44444444444444420e-01 4.58333333333333315e-01 4.72222222222222210e-01 4.86111111111111105e-01 5.00000000000000000e-01 + 1.57079632679489656e+00 1.55691184868293275e+00 1.54303360701525194e+00 1.52915772412750273e+00 1.51528032235533483e+00 1.50139752403439686e+00 1.48750545150033808e+00 1.47360022708880711e+00 1.45967797313545344e+00 1.44573481197592613e+00 1.43176686594587399e+00 1.41777025738094609e+00 1.40374110861679124e+00 1.38967554198905918e+00 1.37556967983339828e+00 1.36141964448545805e+00 1.34722155828088708e+00 1.33297154355533487e+00 1.31866572264445048e+00 1.30430021788388251e+00 1.28987115160928023e+00 1.27537464615629270e+00 1.26080682386056875e+00 1.24616380705775787e+00 1.23144171808350888e+00 1.21663667927347063e+00 1.20174481296329239e+00 1.18676224148862319e+00 1.17168508718511211e+00 1.15650947238840818e+00 1.14123151943416024e+00 1.12584735065801733e+00 1.11035308839562874e+00 1.09474485498264351e+00 1.07901877275471048e+00 1.06317096404747891e+00 1.04719755119659785e+00 +-1.00000000000000000e+00 -9.99411380094997748e-01 -9.99101952023264350e-01 -9.99071715784799697e-01 -9.99320671379603787e-01 -9.99848818807676620e-01 -1.00065615806901831e+00 -1.00174268916362874e+00 -1.00310841209150792e+00 -1.00475332685265606e+00 -1.00667743344707272e+00 -1.00888073187475835e+00 -1.01136322213571272e+00 -1.01412490422993584e+00 -1.01716577815742770e+00 -1.02048584391818853e+00 -1.02408510151221788e+00 -1.02796355093951619e+00 -1.03212119220008325e+00 -1.03655802529391905e+00 -1.04127405022102359e+00 -1.04626926698139711e+00 -1.05154367557503914e+00 -1.05709727600195014e+00 -1.06293006826212988e+00 -1.06904205235557836e+00 -1.07543322828229559e+00 -1.08210359604228179e+00 -1.08905315563553651e+00 -1.09628190706206019e+00 -1.10378985032185262e+00 -1.11157698541491379e+00 -1.11964331234124370e+00 -1.12798883110084258e+00 -1.13661354169370998e+00 -1.14551744411984635e+00 -1.15470053837925146e+00 + 5.24315391578369372e-02 3.23297271624847227e-02 1.22279151671325081e-02 -7.87389682821970638e-03 -2.79757088235719209e-02 -4.80775208189241354e-02 -6.81793328142763499e-02 -8.82811448096285367e-02 -1.08382956804980779e-01 -1.28484768800332994e-01 -1.48586580795685208e-01 -1.68688392791037367e-01 -1.88790204786389637e-01 -2.08892016781741852e-01 -2.28993828777094011e-01 -2.49095640772446225e-01 -2.69197452767798495e-01 -2.89299264763150710e-01 -3.09401076758502924e-01 -3.29502888753855139e-01 -3.49604700749207353e-01 -3.69706512744559457e-01 -3.89808324739911671e-01 -4.09910136735263997e-01 -4.30011948730616211e-01 -4.50113760725968426e-01 -4.70215572721320640e-01 -4.90317384716672855e-01 -5.10419196712024958e-01 -5.30521008707377173e-01 -5.50622820702729387e-01 -5.70724632698081602e-01 -5.90826444693433928e-01 -6.10928256688786142e-01 -6.31030068684138357e-01 -6.51131880679490571e-01 -6.71233692674842786e-01 diff --git a/source/module_base/test/data/two_points_periodic.dat b/source/module_base/test/data/two_points_periodic.dat new file mode 100644 index 0000000000..3b9a3d6d4d --- /dev/null +++ b/source/module_base/test/data/two_points_periodic.dat @@ -0,0 +1,7 @@ +periodic periodic + 0.00000000000000000e+00 1.00000000000000000e+00 + 1.00000000000000000e+00 1.00000000000000000e+00 + 0.00000000000000000e+00 2.77777777777777762e-02 5.55555555555555525e-02 8.33333333333333287e-02 1.11111111111111105e-01 1.38888888888888895e-01 1.66666666666666657e-01 1.94444444444444420e-01 2.22222222222222210e-01 2.50000000000000000e-01 2.77777777777777790e-01 3.05555555555555525e-01 3.33333333333333315e-01 3.61111111111111105e-01 3.88888888888888840e-01 4.16666666666666630e-01 4.44444444444444420e-01 4.72222222222222210e-01 5.00000000000000000e-01 5.27777777777777790e-01 5.55555555555555580e-01 5.83333333333333259e-01 6.11111111111111049e-01 6.38888888888888840e-01 6.66666666666666630e-01 6.94444444444444420e-01 7.22222222222222210e-01 7.50000000000000000e-01 7.77777777777777679e-01 8.05555555555555469e-01 8.33333333333333259e-01 8.61111111111111049e-01 8.88888888888888840e-01 9.16666666666666630e-01 9.44444444444444420e-01 9.72222222222222210e-01 1.00000000000000000e+00 + 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 1.00000000000000000e+00 + 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 + 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 diff --git a/source/module_base/test/data/two_points_second_deriv.dat b/source/module_base/test/data/two_points_second_deriv.dat new file mode 100644 index 0000000000..79882a4a18 --- /dev/null +++ b/source/module_base/test/data/two_points_second_deriv.dat @@ -0,0 +1,7 @@ +second_deriv:5.00000000000000000e-01 second_deriv:2.50000000000000000e-01 + 1.00000000000000000e+00 4.00000000000000000e+00 + 1.00000000000000000e+00 2.00000000000000000e+00 + 1.00000000000000000e+00 1.08333333333333326e+00 1.16666666666666674e+00 1.25000000000000000e+00 1.33333333333333326e+00 1.41666666666666652e+00 1.50000000000000000e+00 1.58333333333333326e+00 1.66666666666666652e+00 1.75000000000000000e+00 1.83333333333333326e+00 1.91666666666666652e+00 2.00000000000000000e+00 2.08333333333333304e+00 2.16666666666666652e+00 2.25000000000000000e+00 2.33333333333333304e+00 2.41666666666666652e+00 2.50000000000000000e+00 2.58333333333333304e+00 2.66666666666666652e+00 2.75000000000000000e+00 2.83333333333333304e+00 2.91666666666666652e+00 3.00000000000000000e+00 3.08333333333333304e+00 3.16666666666666652e+00 3.25000000000000000e+00 3.33333333333333304e+00 3.41666666666666652e+00 3.50000000000000000e+00 3.58333333333333304e+00 3.66666666666666652e+00 3.75000000000000000e+00 3.83333333333333304e+00 3.91666666666666652e+00 4.00000000000000000e+00 + 1.00000000000000000e+00 9.77422518004115282e-01 9.58269032921810648e-01 9.42491319444444531e-01 9.30041152263374471e-01 9.20870306069958899e-01 9.14930555555555469e-01 9.12173675411522611e-01 9.12551440329218089e-01 9.16015625000000000e-01 9.22518004115226220e-01 9.32010352366254957e-01 9.44444444444444420e-01 9.59772055041152150e-01 9.77944958847736689e-01 9.98914930555555469e-01 1.02263374485596703e+00 1.04905317644032903e+00 1.07812500000000000e+00 1.10980099022633705e+00 1.14403292181069949e+00 1.18077256944444420e+00 1.21997170781892983e+00 1.26158211162551437e+00 1.30555555555555558e+00 1.35184381430041123e+00 1.40039866255144019e+00 1.45117187500000000e+00 1.50411522633744821e+00 1.55918049125514369e+00 1.61631944444444442e+00 1.67548386059670751e+00 1.73662551440329205e+00 1.79969618055555580e+00 1.86464763374485565e+00 1.93143164866255113e+00 2.00000000000000000e+00 +-2.91666666666666685e-01 -2.50289351851851916e-01 -2.09490740740740727e-01 -1.69270833333333343e-01 -1.29629629629629706e-01 -9.05671296296297335e-02 -5.20833333333333703e-02 -1.41782407407408349e-02 2.31481481481480254e-02 5.98958333333332801e-02 9.60648148148147418e-02 1.31655092592592504e-01 1.66666666666666630e-01 2.01099537037036841e-01 2.34953703703703554e-01 2.68229166666666630e-01 3.00925925925925764e-01 3.33043981481481399e-01 3.64583333333333259e-01 3.95543981481481344e-01 4.25925925925925930e-01 4.55729166666666741e-01 4.84953703703703665e-01 5.13599537037036979e-01 5.41666666666666741e-01 5.69155092592592338e-01 5.96064814814814659e-01 6.22395833333333259e-01 6.48148148148147918e-01 6.73321759259259189e-01 6.97916666666666630e-01 7.21932870370370239e-01 7.45370370370370239e-01 7.68229166666666630e-01 7.90509259259259189e-01 8.12210648148148140e-01 8.33333333333333370e-01 + 4.99999999999999944e-01 4.93055555555555525e-01 4.86111111111111049e-01 4.79166666666666630e-01 4.72222222222222210e-01 4.65277777777777790e-01 4.58333333333333315e-01 4.51388888888888895e-01 4.44444444444444475e-01 4.37500000000000000e-01 4.30555555555555580e-01 4.23611111111111160e-01 4.16666666666666685e-01 4.09722222222222265e-01 4.02777777777777846e-01 3.95833333333333370e-01 3.88888888888888951e-01 3.81944444444444531e-01 3.75000000000000056e-01 3.68055555555555636e-01 3.61111111111111160e-01 3.54166666666666741e-01 3.47222222222222321e-01 3.40277777777777901e-01 3.33333333333333426e-01 3.26388888888889006e-01 3.19444444444444531e-01 3.12500000000000111e-01 3.05555555555555691e-01 2.98611111111111271e-01 2.91666666666666796e-01 2.84722222222222376e-01 2.77777777777777957e-01 2.70833333333333481e-01 2.63888888888889062e-01 2.56944444444444642e-01 2.50000000000000167e-01 diff --git a/source/module_basis/module_nao/numerical_radial.cpp b/source/module_basis/module_nao/numerical_radial.cpp index 5573fa0ddf..703dece57c 100644 --- a/source/module_basis/module_nao/numerical_radial.cpp +++ b/source/module_basis/module_nao/numerical_radial.cpp @@ -237,8 +237,7 @@ void NumericalRadial::set_grid(const bool for_r_space, const int ngrid, const do #endif // cubic spline interpolation - ModuleBase::CubicSpline cubspl; - cubspl.build(ngrid_tbu, grid_tbu, value_tbu); // not-a-knot boundary condition + ModuleBase::CubicSpline cubspl(ngrid_tbu, grid_tbu, value_tbu); // not-a-knot boundary condition double* grid_new = new double[ngrid]; double* value_new = new double[ngrid]; diff --git a/source/module_basis/module_nao/projgen.cpp b/source/module_basis/module_nao/projgen.cpp index 568551ce5e..dc8d880019 100644 --- a/source/module_basis/module_nao/projgen.cpp +++ b/source/module_basis/module_nao/projgen.cpp @@ -116,16 +116,14 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc }; // cubic spline interpolation - ModuleBase::CubicSpline cubspl; - cubspl.build(nr_proj, r, chi); + ModuleBase::CubicSpline cubspl(nr_proj, r, chi); std::vector dchi(nr_proj); cubspl.eval(nr_proj, r, nullptr, dchi.data()); // function for calculating the overlap between dalpha and dchi auto overlap_dalpha_dchi = [&]() { // calculate dalpha first - ModuleBase::CubicSpline cubspl_alpha; - cubspl_alpha.build(nr_proj, r, alpha.data()); + ModuleBase::CubicSpline cubspl_alpha(nr_proj, r, alpha.data()); std::vector dalpha(nr_proj); cubspl_alpha.eval(nr_proj, r, nullptr, dalpha.data()); for (int i = 0; i < nr_proj; i++) diff --git a/source/module_basis/module_nao/two_center_table.cpp b/source/module_basis/module_nao/two_center_table.cpp index ada89578d5..969e118699 100644 --- a/source/module_basis/module_nao/two_center_table.cpp +++ b/source/module_basis/module_nao/two_center_table.cpp @@ -222,8 +222,10 @@ void TwoCenterTable::_tabulate(const NumericalRadial* it1, const NumericalRadial // than two spherical Bessel transforms. By doing so, we achieve a good // consistency between the table and its derivative during interpolation. using ModuleBase::CubicSpline; - CubicSpline::build(nr_, rgrid_, table_.inner_most_ptr(itab), dtable_.inner_most_ptr(itab), - CubicSpline::BoundaryCondition::first_deriv, CubicSpline::BoundaryCondition::first_deriv, - 0.0, 0.0); + CubicSpline::build( + nr_, rgrid_, table_.inner_most_ptr(itab), + {CubicSpline::BoundaryType::first_deriv, 0.0}, + {CubicSpline::BoundaryType::first_deriv, 0.0}, + dtable_.inner_most_ptr(itab)); }