Add simple library for automatic differentiation

No build system or automatic test cases at this point.

AutoDiff.hpp

@@ -0,0 +1,310 @@
+/*===========================================================================
+//
+// File: AutoDiff.hpp
+//
+// Created: 2013-04-29 10:51:23+0200
+//
+// Authors: Knut-Andreas Lie      <Knut-Andreas.Lie@sintef.no>
+//          Halvor M. Nilsen      <HalvorMoll.Nilsen@sintef.no>
+//          Atgeirr F. Rasmussen  <atgeirr@sintef.no>
+//          Xavier Raynaud        <Xavier.Raynaud@sintef.no>
+//          Bård Skaflestad       <Bard.Skaflestad@sintef.no>
+//
+//==========================================================================*/
+
+
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <cmath>
+
+#ifndef OPM_AUTODIFF_HPP_HEADER
+#define OPM_AUTODIFF_HPP_HEADER
+
+namespace AutoDiff {
+    template <typename Scalar>
+    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 <class Ostream>
+        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 <class Ostream, typename Scalar>
+    Ostream&
+    operator<<(Ostream& os, const Forward<Scalar>& fw)
+    {
+        return fw.print(os);
+    }
+
+    template <typename Scalar>
+    Forward<Scalar>
+    operator +(const Forward<Scalar>& lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Forward<Scalar> ret = lhs;
+        ret += rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator +(const T                lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Forward<Scalar> ret = rhs;
+        ret += Scalar(lhs);
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator +(const Forward<Scalar>& lhs,
+               const T                rhs)
+    {
+        Forward<Scalar> ret = lhs;
+        ret += Scalar(rhs);
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    Forward<Scalar>
+    operator -(const Forward<Scalar>& lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Forward<Scalar> ret = lhs;
+        ret -= rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator -(const T                lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Forward<Scalar> ret(Scalar(lhs) - rhs.val(), -rhs.der());
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator -(const Forward<Scalar>& lhs,
+               const T                rhs)
+    {
+        Forward<Scalar> ret = lhs;
+        ret -= Scalar(rhs);
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    Forward<Scalar>
+    operator *(const Forward<Scalar>& lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Forward<Scalar> ret = lhs;
+        ret *= rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator *(const T                lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Forward<Scalar> ret = rhs;
+        ret *= Scalar(lhs);
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator *(const Forward<Scalar>& lhs,
+               const T                rhs)
+    {
+        Forward<Scalar> ret = lhs;
+        ret *= Scalar(rhs);
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    Forward<Scalar>
+    operator /(const Forward<Scalar>& lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Forward<Scalar> ret = lhs;
+        ret /= rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator /(const T                lhs,
+               const Forward<Scalar>& rhs)
+    {
+        Scalar a =  Scalar(lhs) / rhs.val();
+        Scalar b = -Scalar(lhs) / (rhs.val() * rhs.val());
+
+        Forward<Scalar> ret(a, b);
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    Forward<Scalar>
+    operator /(const Forward<Scalar>& lhs,
+               const T                rhs)
+    {
+        Scalar a = rhs.val() / Scalar(lhs);
+        Scalar b = rhs.der() / Scalar(lhs);
+
+        Forward<Scalar> ret(a, b);
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    Forward<Scalar>
+    cos(const Forward<Scalar>& x)
+    {
+        Forward<Scalar> ret( std::cos(x.val()),
+                            -std::sin(x.val()) * x.der());
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    Forward<Scalar>
+    sqrt(const Forward<Scalar>& x)
+    {
+        Scalar a = std::sqrt(x.val());
+        Scalar b = Scalar(1.0) / (Scalar(2.0) * a);
+        Forward<Scalar> ret(a, b * x.der());
+
+        return ret;
+    }
+} // namespace AutoDiff
+
+namespace std {
+    using AutoDiff::cos;
+    using AutoDiff::sqrt;
+}
+
+#endif  /* OPM_AUTODIFF_HPP_HEADER */

find_zero.cpp

@@ -0,0 +1,101 @@
+/*===========================================================================
+//
+// File: find_zero.cpp
+//
+// Created: 2013-04-29 11:58:29+0200
+//
+// Authors: Knut-Andreas Lie      <Knut-Andreas.Lie@sintef.no>
+//          Halvor M. Nilsen      <HalvorMoll.Nilsen@sintef.no>
+//          Atgeirr F. Rasmussen  <atgeirr@sintef.no>
+//          Xavier Raynaud        <Xavier.Raynaud@sintef.no>
+//          Bård Skaflestad       <Bard.Skaflestad@sintef.no>
+//
+//==========================================================================*/
+
+
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "AutoDiff.hpp"
+
+#include <iostream>
+#include <cmath>
+
+struct Func
+{
+    template <typename T>
+    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 ErrorPolicy = ThrowOnError>
+class Newton
+{
+public:
+    /// Implements a scalar Newton solve.
+    template <class Functor>
+    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<double> xfad(x);
+            AutoDiff::Forward<double> 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';
+}

test_ad.cpp

@@ -0,0 +1,67 @@
+/*===========================================================================
+//
+// File: test_ad.cpp
+//
+// Created: 2013-04-29 11:12:34+0200
+//
+// Authors: Knut-Andreas Lie      <Knut-Andreas.Lie@sintef.no>
+//          Halvor M. Nilsen      <HalvorMoll.Nilsen@sintef.no>
+//          Atgeirr F. Rasmussen  <atgeirr@sintef.no>
+//          Xavier Raynaud        <Xavier.Raynaud@sintef.no>
+//          Bård Skaflestad       <Bard.Skaflestad@sintef.no>
+//
+//==========================================================================*/
+
+
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "AutoDiff.hpp"
+
+#include <iostream>
+
+int
+main()
+{
+    AutoDiff::Forward<double> a(0.0);
+    AutoDiff::Forward<double> 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<double>(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<double>(1);
+    std::cout << "a:     " << a << '\n';
+    std::cout << "a - 1: " << (a - 1) << '\n';
+    std::cout << "a - b: " << (a - b) << '\n';
+}