diff --git a/AutoDiff.hpp b/AutoDiff.hpp new file mode 100644 index 000000000..e7df97634 --- /dev/null +++ b/AutoDiff.hpp @@ -0,0 +1,310 @@ +/*=========================================================================== +// +// File: AutoDiff.hpp +// +// Created: 2013-04-29 10:51:23+0200 +// +// Authors: Knut-Andreas Lie +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Xavier Raynaud +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2013 Statoil ASA. + + 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 + +#ifndef OPM_AUTODIFF_HPP_HEADER +#define OPM_AUTODIFF_HPP_HEADER + +namespace AutoDiff { + template + class Forward { + public: + explicit Forward(const Scalar& x) + : x_(x), dx_(Scalar(1)) + {} + + Forward(const Scalar x, const Scalar dx) + : x_(x), dx_(dx) + {} + + Forward& + operator +=(const Scalar& rhs) + { + x_ += rhs; + + return *this; + } + + Forward& + operator +=(const Forward& rhs) + { + x_ += rhs.x_; + dx_ += rhs.dx_; + + return *this; + } + + Forward& + operator -=(const Scalar& rhs) + { + x_ -= rhs; + + return *this; + } + + Forward& + operator -=(const Forward& rhs) + { + x_ -= rhs.x_; + dx_ -= rhs.dx_; + + return *this; + } + + Forward& + operator *=(const Scalar& rhs) + { + x_ *= rhs; + dx_ *= rhs; + + return *this; + } + + Forward& + operator *=(const Forward& rhs) + { + x_ *= rhs.x_; + dx_ *= dx_*rhs.x_ + x_*rhs.dx_; + + return *this; + } + + Forward& + operator /=(const Scalar& rhs) + { + x_ /= rhs; + dx_ /= rhs; + + return *this; + } + + Forward& + operator /=(const Forward& rhs) + { + x_ /= rhs.x_; + dx_ = (dx_*rhs.x_ - x_*rhs.dx_) / (rhs.x_ * rhs.x_); + + return *this; + } + + template + Ostream& + print(Ostream& os) const + { + os << "(x,dx) = (" << x_ << ',' << dx_ << ")"; + + return os; + } + + const Scalar val() const { return x_ ; } + const Scalar der() const { return dx_; } + + private: + Scalar x_ ; + Scalar dx_; + + }; + + template + Ostream& + operator<<(Ostream& os, const Forward& fw) + { + return fw.print(os); + } + + template + Forward + operator +(const Forward& lhs, + const Forward& rhs) + { + Forward ret = lhs; + ret += rhs; + + return ret; + } + + template + Forward + operator +(const T lhs, + const Forward& rhs) + { + Forward ret = rhs; + ret += Scalar(lhs); + + return ret; + } + + template + Forward + operator +(const Forward& lhs, + const T rhs) + { + Forward ret = lhs; + ret += Scalar(rhs); + + return ret; + } + + template + Forward + operator -(const Forward& lhs, + const Forward& rhs) + { + Forward ret = lhs; + ret -= rhs; + + return ret; + } + + template + Forward + operator -(const T lhs, + const Forward& rhs) + { + Forward ret(Scalar(lhs) - rhs.val(), -rhs.der()); + + return ret; + } + + template + Forward + operator -(const Forward& lhs, + const T rhs) + { + Forward ret = lhs; + ret -= Scalar(rhs); + + return ret; + } + + template + Forward + operator *(const Forward& lhs, + const Forward& rhs) + { + Forward ret = lhs; + ret *= rhs; + + return ret; + } + + template + Forward + operator *(const T lhs, + const Forward& rhs) + { + Forward ret = rhs; + ret *= Scalar(lhs); + + return ret; + } + + template + Forward + operator *(const Forward& lhs, + const T rhs) + { + Forward ret = lhs; + ret *= Scalar(rhs); + + return ret; + } + + template + Forward + operator /(const Forward& lhs, + const Forward& rhs) + { + Forward ret = lhs; + ret /= rhs; + + return ret; + } + + template + Forward + operator /(const T lhs, + const Forward& rhs) + { + Scalar a = Scalar(lhs) / rhs.val(); + Scalar b = -Scalar(lhs) / (rhs.val() * rhs.val()); + + Forward ret(a, b); + + return ret; + } + + template + Forward + operator /(const Forward& lhs, + const T rhs) + { + Scalar a = rhs.val() / Scalar(lhs); + Scalar b = rhs.der() / Scalar(lhs); + + Forward ret(a, b); + + return ret; + } + + template + Forward + cos(const Forward& x) + { + Forward ret( std::cos(x.val()), + -std::sin(x.val()) * x.der()); + + return ret; + } + + template + Forward + sqrt(const Forward& x) + { + Scalar a = std::sqrt(x.val()); + Scalar b = Scalar(1.0) / (Scalar(2.0) * a); + Forward ret(a, b * x.der()); + + return ret; + } +} // namespace AutoDiff + +namespace std { + using AutoDiff::cos; + using AutoDiff::sqrt; +} + +#endif /* OPM_AUTODIFF_HPP_HEADER */ diff --git a/find_zero.cpp b/find_zero.cpp new file mode 100644 index 000000000..fea87f60b --- /dev/null +++ b/find_zero.cpp @@ -0,0 +1,101 @@ +/*=========================================================================== +// +// File: find_zero.cpp +// +// Created: 2013-04-29 11:58:29+0200 +// +// Authors: Knut-Andreas Lie +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Xavier Raynaud +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2013 Statoil ASA. + + 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 "AutoDiff.hpp" + +#include +#include + +struct Func +{ + template + T operator()(T x) const + { +#if 1 + T r = std::sqrt(std::cos(x * x) + x) - 1.2; + return r; +#else + return x; + // const int n = 6; + // double xv[6] = { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 }; + // double yv[6] = { -0.5, -0.3, -0.1, 0.1, 0.3, 0.5 }; + // int interv = -1; + // for (int i = 0; i < n; ++i) { + // if (x < xv[i]) { + // interv = i - 1; + // break; + // } + // } + // T t = (x - xv[interv])/(xv[interv+1] - xv[interv]); + // return (1.0 - t)*yv[interv] + t*yv[interv+1]; +#endif + } +}; + +// template +class Newton +{ +public: + /// Implements a scalar Newton solve. + template + inline static double solve(const Functor& f, + const double initial_guess, + const int max_iter, + const double tolerance, + int& iterations_used) + { + double x = initial_guess; + iterations_used = 0; + while (std::abs(f(x)) > tolerance && ++iterations_used < max_iter) { + AutoDiff::Forward xfad(x); + AutoDiff::Forward rfad = f(xfad); + x = x - rfad.val()/rfad.der(); + } + return x; + } +}; + + +int main() +{ + int iter = 0; + const double atol = 1.0e-13; + const double soln = Newton::solve(Func(), 0.1, 30, atol, iter); + + std::cout.precision(16); + std::cout << "Solution is: " << soln + << " using " << iter << " iterations." << '\n'; + std::cout << " f(x) = " << Func()(soln) << '\n'; +} diff --git a/test_ad.cpp b/test_ad.cpp new file mode 100644 index 000000000..4b7a6f575 --- /dev/null +++ b/test_ad.cpp @@ -0,0 +1,67 @@ +/*=========================================================================== +// +// File: test_ad.cpp +// +// Created: 2013-04-29 11:12:34+0200 +// +// Authors: Knut-Andreas Lie +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Xavier Raynaud +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2013 Statoil ASA. + + 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 "AutoDiff.hpp" + +#include + +int +main() +{ + AutoDiff::Forward a(0.0); + AutoDiff::Forward b(1.0); + + std::cout << "a: " << a << '\n'; + std::cout << "b: " << b << '\n'; + std::cout << "a + b: " << a + b << '\n'; + + a = b; + std::cout << "a: " << a << '\n'; + std::cout << "b: " << b << '\n'; + std::cout << "a + b: " << a + b << '\n'; + + a = AutoDiff::Forward(0.0); + std::cout << "a: " << a << '\n'; + + a += 1; + std::cout << "a: " << a << '\n'; + std::cout << "a + 1: " << (a + 1) << '\n'; + std::cout << "1 + a: " << (1.0f + a) << '\n'; + + a = AutoDiff::Forward(1); + std::cout << "a: " << a << '\n'; + std::cout << "a - 1: " << (a - 1) << '\n'; + std::cout << "a - b: " << (a - b) << '\n'; +}