From 2dc1ec345111034986eed347aac6943dd740e96f Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Tue, 8 Aug 2023 23:38:48 +0800 Subject: [PATCH] new CubicSpline interface (#2800) Co-authored-by: Wenfei Li <38569667+wenfei-li@users.noreply.github.com> --- source/module_base/cubic_spline.cpp | 217 ++++++++++++------ source/module_base/cubic_spline.h | 199 +++++++++++----- source/module_base/test/cubic_spline_test.cpp | 33 ++- .../module_nao/numerical_radial.cpp | 2 +- 4 files changed, 314 insertions(+), 137 deletions(-) diff --git a/source/module_base/cubic_spline.cpp b/source/module_base/cubic_spline.cpp index a6ea536211..230ae56d74 100644 --- a/source/module_base/cubic_spline.cpp +++ b/source/module_base/cubic_spline.cpp @@ -9,9 +9,40 @@ namespace ModuleBase { -CubicSpline::~CubicSpline() +CubicSpline::CubicSpline(CubicSpline const& other) { - 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)); +} + +CubicSpline& CubicSpline::operator=(CubicSpline const& other) +{ + if (this != &other) + { + 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)); + } + + return *this; } void CubicSpline::cleanup() @@ -25,13 +56,12 @@ void CubicSpline::cleanup() s_ = nullptr; } -void CubicSpline::sanity_check(const int n, - const double* const x, - const double* const y, - BoundaryCondition bc_start, - BoundaryCondition bc_end) +void CubicSpline::check_build(const int n, + const double* const x, + const double* const y, + BoundaryCondition bc_start, + BoundaryCondition bc_end) { - assert(n > 1); // periodic boundary condition must apply to both ends @@ -49,30 +79,38 @@ void CubicSpline::sanity_check(const int n, } } +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) +{ + 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 + + // check that x_interp is in the range of the interpolant + assert(std::all_of(x_interp, x_interp + n_interp, + [n, x](const double xi) { return xi >= x[0] && xi <= x[n - 1]; })); +} + 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) { - sanity_check(n, x, y, bc_start, bc_end); - - cleanup(); - - n_ = n; - x_ = new double[n]; - y_ = new double[n]; - std::memcpy(x_, x, sizeof(double) * n); - std::memcpy(y_, y, sizeof(double) * n); - - // to be computed - s_ = new double[n]; + check_build(n, x, 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; + s[0] = s[1] = 0.0; } 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 @@ -80,9 +118,9 @@ void CubicSpline::build(const int n, 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]); + 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 { @@ -94,18 +132,10 @@ void CubicSpline::build(const int n, double* subdiag = new double[n - 1]; double* supdiag = new double[n - 1]; - is_uniform_ = true; - double dx_avg = (x[n - 1] - x[0]) / (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]; - - if (std::abs(dx[i] - dx_avg) > uniform_thr_) - { - is_uniform_ = false; - } } // common part of the tridiagonal linear system @@ -114,20 +144,18 @@ void CubicSpline::build(const int n, 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]); + s[i] = 3.0 * (dd[i - 1] * dx[i] + dd[i] * dx[i - 1]); } if (bc_start == BoundaryCondition::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, supdiag, subdiag, s_); - s_[n - 1] = s_[0]; + s[0] = 3.0 * (dd[0] * dx[n - 2] + dd[n - 2] * dx[0]); + solve_cyctri(n - 1, diag, supdiag, subdiag, s); + s[n - 1] = s[0]; } else { @@ -136,17 +164,17 @@ void CubicSpline::build(const int n, case BoundaryCondition::first_deriv: diag[0] = 1.0 * dx[0]; supdiag[0] = 0.0; - s_[0] = deriv_start * dx[0]; + s[0] = deriv_start * 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]; + s[0] = (3.0 * dd[0] - 0.5 * deriv_start * 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]); + s[0] = (dd[0] * dx[1] * (dx[0] + 2 * (x[2] - x[0])) + dd[1] * dx[0] * dx[0]) / (x[2] - x[0]); } switch (bc_end) @@ -154,17 +182,17 @@ void CubicSpline::build(const int n, 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]; + s[n - 1] = deriv_end * 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]; + s[n - 1] = (3.0 * dd[n - 2] + 0.5 * deriv_end * 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])) + 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]); } @@ -174,7 +202,7 @@ void CubicSpline::build(const int n, int INFO = 0; int N = n; - dgtsv_(&N, &NRHS, subdiag, diag, supdiag, s_, &LDB, &INFO); + dgtsv_(&N, &NRHS, subdiag, diag, supdiag, s, &LDB, &INFO); } delete[] diag; @@ -185,41 +213,98 @@ void CubicSpline::build(const int n, } } -void CubicSpline::interp(const int n, const double* const x, double* const y, double* const dy) +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) { - assert(x_ && y_ && s_); // make sure the interpolant exists - assert(y || dy); // make sure at least one of y or dy is not null + cleanup(); - // check that x is in the range of the interpolant - assert(std::all_of(x, x + n, [this](double x) -> bool { return x >= x_[0] && x <= x_[n_ - 1]; })); + n_ = n; + x_ = new double[n]; + y_ = new double[n]; + s_ = new double[n]; - std::function search; - if (is_uniform_) + 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) +{ + if (is_uniform != 0 && is_uniform != 1) { - double dx = x_[1] - x_[0]; - search = [this, dx](double x) -> int { return x == x_[n_ - 1] ? n_ - 2 : x / dx; }; + 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; }); + } + + if (is_uniform) + { + double dx = x[1] - x[0]; + return [dx, n, x](double xi) -> int { return xi == x[n - 1] ? n - 2 : xi / dx; }; } else { - search = [this](double x) -> int { return (std::upper_bound(x_, x_ + n_, x) - x_) - 1; }; + return [n, x](double xi) -> int { return (std::upper_bound(x, x + n, xi) - x) - 1; }; } +} - void (*eval)(double, double, double, double, double, double*, double*) = y ? (dy ? _eval_y_dy : _eval_y) : _eval_dy; +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; ++i) + for (int i = 0; i != n_interp; ++i) { - int p = search(x[i]); - double w = x[i] - x_[p]; + 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 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; + 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; - eval(w, c0, c1, c2, c3, y + i, dy + i); + poly_eval(w, c0, c1, c2, c3, y_interp + i, dy_interp + i); } } diff --git a/source/module_base/cubic_spline.h b/source/module_base/cubic_spline.h index 44fec249e2..c025dd7b6e 100644 --- a/source/module_base/cubic_spline.h +++ b/source/module_base/cubic_spline.h @@ -6,13 +6,17 @@ namespace ModuleBase { -//! A class that performs cubic spline interplations. /*! - * 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]. + * @brief A class that performs cubic spline interplations. * - * Usage: + * 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]. + * + * 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. + * + * Usage-1: interpolant object * * CubicSpline cubspl; * @@ -36,31 +40,42 @@ namespace ModuleBase * CubicSpline::BoundaryCondition::periodic ); * * // calculate the values of interpolant at x_interp[i] - * cubspl.interp(n_interp, x_interp, y_interp); + * cubspl.eval(n_interp, x_interp, y_interp); * * // calculate the values & derivatives of interpolant at x_interp[i] - * cubspl.interp(n_interp, x_interp, y_interp, dy_interp); - * */ + * 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); + * + * // evaluates the interpolant with knots, values & derivatives + * CubicSpline::eval(n, x, y, s, n_interp, x_interp, y_interp, dy_interp); + * + * */ class CubicSpline { public: CubicSpline(){}; - CubicSpline(CubicSpline const&) = delete; - CubicSpline& operator=(CubicSpline const&) = delete; + CubicSpline(CubicSpline const&); + CubicSpline& operator=(CubicSpline const&); - ~CubicSpline(); + ~CubicSpline() { cleanup(); } - //! Boundary condition 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. - * */ + * @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 { not_a_knot, @@ -69,36 +84,41 @@ class CubicSpline periodic }; - //! Builds the interpolant. /*! - * This function computes all cubic spline polynomial coefficients from given - * data points and boundary conditions. + * @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 the start - BoundaryCondition bc_end = BoundaryCondition::not_a_knot, //!< [in] boundary condition at the end + 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 ); - //! Calculates the values of interpolant. /*! - * 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. + * @brief Evaluates the interpolant. + * + * 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. + * + * If y_interp or dy_interp is nullptr, the corresponding values are + * not calculated. They must not be nullptr at the same time. * - * If y_interp or dy_interp is nullptr, the corresponding values are - * not calculated. y_interp and dy_interp must not be both nullptr. + * @note the interpolant must be built before calling this function. * */ - void interp(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 + 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 ); /// knots of the interpolant @@ -110,6 +130,29 @@ class CubicSpline /// 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 + ); + + 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 + ); + private: //! number of data points int n_ = 0; @@ -134,12 +177,22 @@ class CubicSpline * */ double uniform_thr_ = 1e-15; - //! Checks whether the input arguments of build() are valid - void sanity_check(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 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. /*! @@ -160,32 +213,58 @@ class CubicSpline * * @note D, L, U are all modified in this function, so use with care! * */ - 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. + 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. ); //! Wipes off the interpolant (if any) and deallocates memories. void cleanup(); - static void _eval_y(double w, double c0, double c1, double c2, double c3, double* y, double*) + /// 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) { - *y = ((c3 * w + c2) * w + c1) * w + c0; - } + if (EvalY) + { + *y = ((c3 * w + c2) * w + c1) * w + c0; + } - static void _eval_dy(double w, double, double c1, double c2, double c3, double*, double* dy) - { - *dy = (3.0 * c3 * w + 2.0 * c2) * w + c1; + if (EvalDy) + { + *dy = (3.0 * c3 * w + 2.0 * c2) * w + c1; + } } - static void _eval_y_dy(double w, double c0, double c1, double c2, double c3, double* y, double* dy) - { - *y = ((c3 * w + c2) * w + c1) * w + c0; - *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 + ); }; }; // namespace ModuleBase diff --git a/source/module_base/test/cubic_spline_test.cpp b/source/module_base/test/cubic_spline_test.cpp index 49058a4a83..e175ba7330 100644 --- a/source/module_base/test/cubic_spline_test.cpp +++ b/source/module_base/test/cubic_spline_test.cpp @@ -22,7 +22,7 @@ using ModuleBase::PI; * Constructs the cubic spline interpolant from given * data points and boundary condition. * - * - get + * - eval * Evaluates the interpolant at certain points. * */ class CubicSplineTest : public ::testing::Test @@ -37,6 +37,7 @@ class CubicSplineTest : public ::testing::Test 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_ @@ -50,6 +51,7 @@ 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]; @@ -59,6 +61,7 @@ void CubicSplineTest::TearDown() { delete[] x_; delete[] y_; + delete[] s_; delete[] x_interp_; delete[] y_interp_; delete[] y_ref_; @@ -75,6 +78,7 @@ TEST_F(CubicSplineTest, NotAKnot) y_[i] = std::sin(x_[i]); } + // interpolant object cubspl.build(n, x_, y_, CubicSpline::BoundaryCondition::not_a_knot, CubicSpline::BoundaryCondition::not_a_knot); int ni = 6; @@ -92,7 +96,16 @@ TEST_F(CubicSplineTest, NotAKnot) y_ref_[4] = 0.1510153464180796; y_ref_[5] = 0.1411200080598672; - cubspl.interp(ni, x_interp_, y_interp_); + 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); @@ -125,7 +138,7 @@ TEST_F(CubicSplineTest, PeriodicAndUniform) y_ref_[3] = 1.4356324132349871e-04; y_ref_[4] = 1.0000000000000000e+00; - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); @@ -166,7 +179,7 @@ TEST_F(CubicSplineTest, FirstDeriv) y_ref_[4] = 0.0798871927951471; y_ref_[5] = 0.0497870683678639; - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); @@ -207,7 +220,7 @@ TEST_F(CubicSplineTest, SecondDeriv) y_ref_[4] = 0.0788753653329337; y_ref_[5] = 0.0497870683678639; - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); @@ -245,7 +258,7 @@ TEST_F(CubicSplineTest, TwoPoints) y_ref_[3] = 4.159700000000001; y_ref_[4] = 4.33; - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); @@ -266,7 +279,7 @@ TEST_F(CubicSplineTest, TwoPoints) y_ref_[3] = 4.074050000000001; y_ref_[4] = 4.33; - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); @@ -276,7 +289,7 @@ TEST_F(CubicSplineTest, TwoPoints) y_[1] = y_[0]; cubspl.build(2, x_, y_, CubicSpline::BoundaryCondition::periodic, CubicSpline::BoundaryCondition::periodic); - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], 2.33, tol); @@ -313,7 +326,7 @@ TEST_F(CubicSplineTest, ThreePoints) y_ref_[3] = 4.3875; y_ref_[4] = 4.8; - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); @@ -329,7 +342,7 @@ TEST_F(CubicSplineTest, ThreePoints) y_ref_[3] = 1.2361111111111112; y_ref_[4] = 1.2; - cubspl.interp(ni, x_interp_, y_interp_); + cubspl.eval(ni, x_interp_, y_interp_); for (int i = 0; i != ni; ++i) { EXPECT_NEAR(y_interp_[i], y_ref_[i], tol); diff --git a/source/module_basis/module_nao/numerical_radial.cpp b/source/module_basis/module_nao/numerical_radial.cpp index 8b1f9993ba..f32c8f207f 100644 --- a/source/module_basis/module_nao/numerical_radial.cpp +++ b/source/module_basis/module_nao/numerical_radial.cpp @@ -240,7 +240,7 @@ void NumericalRadial::set_grid(const bool for_r_space, const int ngrid, const do // grid_end is the first grid point that is strictly greater than grid_tbu[ngrid_tbu-1] double* grid_end = std::upper_bound(grid_new, grid_new + ngrid, grid_tbu[ngrid_tbu - 1]); - cubspl.interp(std::distance(grid_start, grid_end), grid_start, value_new + std::distance(grid_new, grid_start)); + cubspl.eval(std::distance(grid_start, grid_end), grid_start, value_new + std::distance(grid_new, grid_start)); delete[] grid_tbu; delete[] value_tbu;