diff --git a/opm/core/linalg/blas_lapack.h b/opm/core/linalg/blas_lapack.h new file mode 100644 index 000000000..be2b956c4 --- /dev/null +++ b/opm/core/linalg/blas_lapack.h @@ -0,0 +1,148 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BLAS_LAPACK_HEADER_INCLUDED +#define OPM_BLAS_LAPACK_HEADER_INCLUDED + +#ifdef __cplusplus +extern "C" { +#endif + +#if defined(MATLAB_MEX_FILE) && MATLAB_MEX_FILE +#include +#undef MAT_SIZE_T +#define MAT_SIZE_T mwSignedIndex +#endif + +#ifndef MAT_SIZE_T +#define MAT_SIZE_T int +#endif + + +/* C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} */ +void dgemm_(const char *transA , const char *transB , + const MAT_SIZE_T* m, const MAT_SIZE_T* n , const MAT_SIZE_T* k , + const double* a1, const double* A , const MAT_SIZE_T* ldA, + const double* B, const MAT_SIZE_T* ldB, + const double* a2, double* C , const MAT_SIZE_T* ldC); + + +/* C <- a1*A*A' + a2*C *or* C <- a1*A'*A + a2*C */ +void dsyrk_(const char *uplo, const char *trans, + const MAT_SIZE_T *n , const MAT_SIZE_T *k , + const double *a1 , const double *A , const MAT_SIZE_T *ldA, + const double *a2 , double *C , const MAT_SIZE_T *ldC); + + +void dgeqrf_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , + double *A , const MAT_SIZE_T *ld , + double *tau , double *work, + const MAT_SIZE_T *lwork, MAT_SIZE_T *info); + + +void dorgqr_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , const MAT_SIZE_T *k , + double *A , const MAT_SIZE_T *ld , const double *tau, + double *work, const MAT_SIZE_T *lwork, MAT_SIZE_T *info); + +/* A <- LU(A) */ +void dgetrf_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , + double *A , const MAT_SIZE_T *ld, + MAT_SIZE_T *ipiv, MAT_SIZE_T *info); + +/* B <- A \ B, when A is LU(A) from dgetrf() */ +void dgetrs_(const char *trans, const MAT_SIZE_T *n, + const MAT_SIZE_T *nrhs , + const double *A , const MAT_SIZE_T *lda, + const MAT_SIZE_T *ipiv , double *B, + const MAT_SIZE_T *ldb , MAT_SIZE_T *info); + +/* B <- A \ B, tridiagonal A with bands DL, D, DU */ +void dgtsv_(const MAT_SIZE_T *n , + const MAT_SIZE_T *nrhs , + double *DL , + double *D , + double *DU , + double *B , + const MAT_SIZE_T *ldb , + MAT_SIZE_T *info); + +/* B <- A \ B, band matrix A stored in AB with kl subdiagonals, ku superdiagonals */ +void dgbsv_(const MAT_SIZE_T *n , + const MAT_SIZE_T *kl , + const MAT_SIZE_T *ku , + const MAT_SIZE_T *nrhs , + double *AB , + const MAT_SIZE_T *ldab , + MAT_SIZE_T *ipiv , + double *B , + const MAT_SIZE_T *ldb , + MAT_SIZE_T *info); + +/* B <- A \ B, general solver */ +void dgesv_(const MAT_SIZE_T *n, + const MAT_SIZE_T *nrhs , + double *A , + const MAT_SIZE_T *lda , + MAT_SIZE_T *piv , + double *B , + const MAT_SIZE_T *ldb , + MAT_SIZE_T *info); + +/* A <- chol(A) */ +void dpotrf_(const char *uplo, const MAT_SIZE_T *n, + double *A , const MAT_SIZE_T *lda, + MAT_SIZE_T *info); + +/* B <- (A \ (A' \ B)), when A=DPOTRF(A_orig) */ +void dpotrs_(const char *uplo, const MAT_SIZE_T *n , const MAT_SIZE_T *nrhs, + double *A , const MAT_SIZE_T *lda, + double *B , const MAT_SIZE_T *ldb, MAT_SIZE_T *info); + +/* A <- chol(A), packed format. */ +void dpptrf_(const char *uplo, const MAT_SIZE_T *n, + double *Ap , MAT_SIZE_T *info); + +/* A <- (A \ (A' \ eye(n)) when A=DPPTRF(A_orig) (packed format). */ +void dpptri_(const char *uplo, const MAT_SIZE_T *n, + double *Ap , MAT_SIZE_T *info); + +/* y <- a1*op(A)*x + a2*y */ +void dgemv_(const char *trans, + const MAT_SIZE_T *m , const MAT_SIZE_T *n, + const double *a1 , const double *A, const MAT_SIZE_T *ldA , + const double *x, const MAT_SIZE_T *incX, + const double *a2 , double *y, const MAT_SIZE_T *incY); + + +/* y <- a*x + y */ +void daxpy_(const MAT_SIZE_T *n, const double *a, + const double *x, const MAT_SIZE_T *incx, + double *y, const MAT_SIZE_T *incy); + +/* s <- x' * y */ +double ddot_(const MAT_SIZE_T *n, const double *x, const MAT_SIZE_T *incx, + const double *y, const MAT_SIZE_T *incy); + + + +#ifdef __cplusplus +} +#endif + +#endif /* OPM_BLAS_LAPACK_HEADER_INCLUDED */ diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp new file mode 100644 index 000000000..dc7f75e2a --- /dev/null +++ b/opm/core/utility/UniformTableLinear.hpp @@ -0,0 +1,268 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED +#define OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED + +#include +#include +#include +#include +#include + +#include + +namespace Opm { + + /// @brief This class uses linear interpolation to compute the value + /// (and its derivative) of a function f sampled at uniform points. + /// @tparam T the range type of the function (should be an algebraic ring type) + template + class UniformTableLinear + { + public: + /// @brief Default constructor. + UniformTableLinear(); + + /// @brief Construct from vector of y-values. + /// @param xmin the x value corresponding to the first y value. + /// @param xmax the x value corresponding to the last y value. + /// @param y_values vector of range values. + UniformTableLinear(double xmin, + double xmax, + const std::vector& y_values); + + /// @brief Construct from array of y-values. + /// @param xmin the x value corresponding to the first y value. + /// @param xmax the x value corresponding to the last y value. + /// @param y_values array of range values. + /// @param num_y_values the number of values in y_values. + UniformTableLinear(double xmin, + double xmax, + const T* y_values, + int num_y_values); + + /// @brief Get the domain. + /// @return the domain as a pair of doubles. + std::pair domain(); + + /// @brief Rescale the domain. + /// @param new_domain the new domain as a pair of doubles. + void rescaleDomain(std::pair new_domain); + + /// @brief Evaluate the value at x. + /// @param x a domain value + /// @return f(x) + double operator()(const double x) const; + + /// @brief Evaluate the derivative at x. + /// @param x a domain value + /// @return f'(x) + double derivative(const double x) const; + + /// @brief Equality operator. + /// @param other another UniformTableLinear. + /// @return true if they are represented exactly alike. + bool operator==(const UniformTableLinear& other) const; + + /// @brief Policies for how to behave when trying to evaluate outside the domain. + enum RangePolicy {Throw = 0, ClosestValue = 1, Extrapolate = 2}; + + /// @brief Sets the behavioural policy for evaluation to the left of the domain. + /// @param rp the policy + void setLeftPolicy(RangePolicy rp); + + /// @brief Sets the behavioural policy for evaluation to the right of the domain. + /// @param rp the policy + void setRightPolicy(RangePolicy rp); + + protected: + double xmin_; + double xmax_; + double xdelta_; + std::vector y_values_; + RangePolicy left_; + RangePolicy right_; + template + friend std::ostream& operator<<(std::ostream& os, const UniformTableLinear& t); + }; + + + // Member implementations. + + template + inline + UniformTableLinear + ::UniformTableLinear() + : left_(ClosestValue), right_(ClosestValue) + { + } + + template + inline + UniformTableLinear + ::UniformTableLinear(double xmin, + double xmax, + const std::vector& y_values) + : xmin_(xmin), xmax_(xmax), y_values_(y_values), + left_(ClosestValue), right_(ClosestValue) + { + assert(xmax > xmin); + assert(y_values.size() > 1); + xdelta_ = (xmax - xmin)/(y_values.size() - 1); + } + + template + inline + UniformTableLinear + ::UniformTableLinear(double xmin, + double xmax, + const T* y_values, + int num_y_values) + : xmin_(xmin), xmax_(xmax), + y_values_(y_values, y_values + num_y_values), + left_(ClosestValue), right_(ClosestValue) + { + assert(xmax > xmin); + assert(y_values_.size() > 1); + xdelta_ = (xmax - xmin)/(y_values_.size() - 1); + } + + template + inline std::pair + UniformTableLinear + ::domain() + { + return std::make_pair(xmin_, xmax_); + } + + template + inline void + UniformTableLinear + ::rescaleDomain(std::pair new_domain) + { + xmin_ = new_domain.first; + xmax_ = new_domain.second; + xdelta_ = (xmax_ - xmin_)/(y_values_.size() - 1); + } + + template + inline double + UniformTableLinear + ::operator()(const double xparam) const + { + // Implements ClosestValue policy. + double x = std::min(xparam, xmax_); + x = std::max(x, xmin_); + + // Lookup is easy since we are uniform in x. + double pos = (x - xmin_)/xdelta_; + double posi = std::floor(pos); + int left = int(posi); + if (left == int(y_values_.size()) - 1) { + // We are at xmax_ + return y_values_.back(); + } + double w = pos - posi; + return (1.0 - w)*y_values_[left] + w*y_values_[left + 1]; + } + + template + inline double + UniformTableLinear + ::derivative(const double xparam) const + { + // Implements derivative consistent + // with ClosestValue policy for function + double value; + if (xparam > xmax_ || xparam < xmin_) { + value = 0.0; + } else { + double x = std::min(xparam, xmax_); + x = std::max(x, xmin_); + // Lookup is easy since we are uniform in x. + double pos = (x - xmin_)/xdelta_; + double posi = std::floor(pos); + int left = int(posi); + if (left == int(y_values_.size()) - 1) { + // We are at xmax_ + --left; + } + value = (y_values_[left + 1] - y_values_[left])/xdelta_; + } + return value; + } + + + template + inline bool + UniformTableLinear + ::operator==(const UniformTableLinear& other) const + { + return xmin_ == other.xmin_ + && xdelta_ == other.xdelta_ + && y_values_ == other.y_values_ + && left_ == other.left_ + && right_ == other.right_; + } + + template + inline void + UniformTableLinear + ::setLeftPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + OPM_THROW(std::runtime_error, "Only ClosestValue RangePolicy implemented."); + } + left_ = rp; + } + + template + inline void + UniformTableLinear + ::setRightPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + OPM_THROW(std::runtime_error, "Only ClosestValue RangePolicy implemented."); + } + right_ = rp; + } + + + template + inline std::ostream& operator<<(std::ostream& os, const UniformTableLinear& t) + { + int n = t.y_values_.size(); + for (int i = 0; i < n; ++i) { + double f = double(i)/double(n - 1); + os << (1.0 - f)*t.xmin_ + f*t.xmax_ + << " " << t.y_values_[i] << '\n'; + } + return os; + } + + namespace utils + { + using Opm::UniformTableLinear; + } + + +} // namespace Opm + +#endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED diff --git a/opm/core/utility/buildUniformMonotoneTable.hpp b/opm/core/utility/buildUniformMonotoneTable.hpp new file mode 100644 index 000000000..3ebad7c99 --- /dev/null +++ b/opm/core/utility/buildUniformMonotoneTable.hpp @@ -0,0 +1,50 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED +#define OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED + +#include +#include + +namespace Opm { + + template + void buildUniformMonotoneTable(const std::vector& xv, + const std::vector& yv, + const int samples, + UniformTableLinear& table) + { + MonotCubicInterpolator interp(xv, yv); + std::vector uniform_yv(samples); + double xmin = xv[0]; + double xmax = xv.back(); + for (int i = 0; i < samples; ++i) { + double w = double(i)/double(samples - 1); + double x = (1.0 - w)*xmin + w*xmax; + uniform_yv[i] = interp(x); + } + table = UniformTableLinear(xmin, xmax, uniform_yv); + } + +} // namespace Opm + + + +#endif // OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED diff --git a/tests/test_cubic.cpp b/tests/test_cubic.cpp new file mode 100644 index 000000000..ecbb3f686 --- /dev/null +++ b/tests/test_cubic.cpp @@ -0,0 +1,73 @@ +//=========================================================================== +// +// File: monotcubicinterpolator_test.cpp +// +// Created: Tue Dec 8 12:25:30 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $Date$ +// +// $Revision$ +// +//=========================================================================== + +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. + Portions Copyright 2013 Uni Research AS. + + This file is part of The Open Reservoir Simulator Project (OpenRS). + + OpenRS is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenRS is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . +*/ + +#include "config.h" + +/* --- Boost.Test boilerplate --- */ +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE CubicTest +#include +#include + +/* --- our own headers --- */ +#include +using namespace Opm; + +BOOST_AUTO_TEST_SUITE () + +BOOST_AUTO_TEST_CASE (cubic) +{ + const int num_v = 3; + double xv[num_v] = {0.0, 1.0, 2.0}; + double fv[num_v] = {10.0, 21.0, 2.0}; + std::vector x(xv, xv + num_v); + std::vector f(fv, fv + num_v); + MonotCubicInterpolator interp(x, f); + BOOST_REQUIRE_CLOSE (interp.evaluate(-1.0), 10., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(0.0), 10., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(0.0001), 10.0011, 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(0.5), 17.375, 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(1.0), 21., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(2.0), 2., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(4.0), 2., 0.00001); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/tests/test_sparsevector.cpp b/tests/test_sparsevector.cpp new file mode 100644 index 000000000..a6f89be09 --- /dev/null +++ b/tests/test_sparsevector.cpp @@ -0,0 +1,109 @@ +//=========================================================================== +// +// File: sparsevector_test.cpp +// +// Created: Mon Jun 29 21:00:53 2009 +// +// Author(s): Atgeirr F Rasmussen +// BÃ¥rd Skaflestad +// +// $Date$ +// +// $Revision$ +// +//=========================================================================== + +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. + + This file is part of The Open Reservoir Simulator Project (OpenRS). + + OpenRS is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenRS is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE SparseVectorTest +#include + +#include + +using namespace Opm; + +BOOST_AUTO_TEST_CASE(construction_and_queries) +{ + const SparseVector sv1; + BOOST_CHECK(sv1.empty()); + BOOST_CHECK_EQUAL(sv1.size(), 0); + BOOST_CHECK_EQUAL(sv1.nonzeroSize(), 0); + + const int size = 100; + const int num_elem = 9; + const int elem[num_elem] = { 9, 8, 7, 6, 5, 4, 3, 2, 1 }; + const int indices[num_elem] = { 1, 2, 3, 5, 8, 13, 21, 34, 55 }; + const SparseVector sv2(size, elem, elem + num_elem, indices, indices + num_elem); + BOOST_CHECK(!sv2.empty()); + BOOST_CHECK_EQUAL(sv2.size(), size); + BOOST_CHECK_EQUAL(sv2.element(0), 0); + BOOST_CHECK_EQUAL(sv2.element(1), 9); + BOOST_CHECK_EQUAL(sv2.element(2), 8); + BOOST_CHECK_EQUAL(sv2.element(3), 7); + BOOST_CHECK_EQUAL(sv2.element(4), 0); + BOOST_CHECK_EQUAL(sv2.element(5), 6); + BOOST_CHECK_EQUAL(sv2.element(55), 1); + BOOST_CHECK_EQUAL(sv2.element(99), 0); + BOOST_CHECK_EQUAL(sv2.nonzeroSize(), num_elem); + for (int i = 0; i < num_elem; ++i) { + BOOST_CHECK_EQUAL(sv2.nonzeroElement(i), elem[i]); + } + const SparseVector sv2_again(size, elem, elem + num_elem, indices, indices + num_elem); + BOOST_CHECK(sv2 == sv2_again); + SparseVector sv2_append(size, elem, elem + num_elem - 1, indices, indices + num_elem - 1); + BOOST_CHECK_EQUAL(sv2_append.nonzeroSize(), num_elem - 1); + sv2_append.addElement(elem[num_elem - 1], indices[num_elem - 1]); + BOOST_CHECK(sv2 == sv2_append); + SparseVector sv2_append2(size); + for (int i = 0; i < num_elem; ++i) { + sv2_append2.addElement(elem[i], indices[i]); + } + BOOST_CHECK(sv2 == sv2_append2); + sv2_append2.clear(); + SparseVector sv_empty; + BOOST_CHECK(sv2_append2 == sv_empty); + + // Tests that only run in debug mode. +#ifndef NDEBUG + // One element too few. + BOOST_CHECK_THROW(const SparseVector sv3(size, elem, elem + num_elem - 1, indices, indices + num_elem), std::exception); + // One element too many. + BOOST_CHECK_THROW(const SparseVector sv4(size, elem, elem + num_elem, indices, indices + num_elem - 1), std::exception); + // Indices out of range. + BOOST_CHECK_THROW(const SparseVector sv5(4, elem, elem + num_elem, indices, indices + num_elem), std::exception); + // Indices not strictly increasing. Cheating by using the elements as indices. + BOOST_CHECK_THROW(const SparseVector sv5(size, elem, elem + num_elem, elem, elem + num_elem), std::exception); + + // Do not ask for out of range indices. + BOOST_CHECK_THROW(sv1.element(0), std::exception); + BOOST_CHECK_THROW(sv2.element(-1), std::exception); + BOOST_CHECK_THROW(sv2.element(sv2.size()), std::exception); + BOOST_CHECK_THROW(sv2.nonzeroElement(sv2.nonzeroSize()), std::exception); +#endif +} + diff --git a/tests/test_uniformtablelinear.cpp b/tests/test_uniformtablelinear.cpp new file mode 100644 index 000000000..ed7d0ea40 --- /dev/null +++ b/tests/test_uniformtablelinear.cpp @@ -0,0 +1,73 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include + +#if defined(HAVE_DYNAMIC_BOOST_TEST) +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + + +#define BOOST_TEST_MODULE UniformTableLinearTests +#include +#include + + + +BOOST_AUTO_TEST_CASE(table_operations) +{ + // Make a simple table. + double yva[] = { 1.0, -1.0, 3.0, 4.0, 2.0 }; + const int numvals = sizeof(yva)/sizeof(yva[0]); + std::vector yv(yva, yva + numvals); + const double xmin = 1.0; + const double xdelta = 2.5; + const double xmax = xmin + (numvals - 1)*xdelta; + Opm::utils::UniformTableLinear t1(xmin, xmax, yv); + Opm::utils::UniformTableLinear t1_copy1(1.0, 11.0, yv); + Opm::utils::UniformTableLinear t1_copy2(t1); + + // Check equality. + BOOST_CHECK(t1 == t1_copy1); + BOOST_CHECK(t1 == t1_copy2); + + // Check some evaluations. + for (int i = 0; i < numvals; ++i) { + BOOST_CHECK_EQUAL(t1(xmin + i*xdelta), yv[i]); + } + BOOST_CHECK_EQUAL(t1(2.25), 0.0); + BOOST_CHECK_EQUAL(t1(9.75), 3.0); + BOOST_CHECK_EQUAL(t1.derivative(9.75), -2.0/xdelta); + // Until we implement anything but the ClosestValue end policy, we only test that. + BOOST_CHECK_EQUAL(t1(xmin - 1.0), yv[0]); + BOOST_CHECK_EQUAL(t1(xmax + 1.0), yv.back()); + + // Domains. + BOOST_CHECK_EQUAL(t1.domain().first, xmin); + BOOST_CHECK_EQUAL(t1.domain().second, xmin + (numvals-1)*xdelta); + std::pair new_domain(-100.0, 20.0); + t1.rescaleDomain(new_domain); + BOOST_CHECK_EQUAL(t1.domain().first, new_domain.first); + BOOST_CHECK_EQUAL(t1.domain().second, new_domain.second); + for (int i = 0; i < numvals; ++i) { + BOOST_CHECK_EQUAL(t1(-100.0 + i*120.0/(double(numvals - 1))), yv[i]); + } + BOOST_CHECK_EQUAL(t1(-85.0), 0.0); + BOOST_CHECK(std::fabs(t1.derivative(0.0) + 2.0/30.0) < 1e-14); +}