Initial work in progress on block AD class.

AutoDiffBlock.hpp

@@ -0,0 +1,157 @@
+/*
+  Copyright 2013 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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
+#define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
+
+#include "AutoDiff.hpp"
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+#include <vector>
+#include <cassert>
+
+namespace AutoDiff
+{
+
+    template <typename Scalar>
+    class ForwardBlock
+    {
+    public:
+        /// Underlying types for scalar vectors and jacobians.
+        typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
+        typedef Eigen::SparseMatrix<Scalar> M;
+
+        /// Named constructor pattern used here.
+        static ForwardBlock constant(const int index, const V& val, const std::vector<int>& blocksizes)
+        {
+            std::vector<M> jac;
+            const int num_elem = val.size();
+            const int num_blocks = blocksizes.size();
+            // For constants, all jacobian blocks are zero.
+            jac.resize(num_blocks);
+            for (int i = 0; i < num_blocks; ++i) {
+                jac[i] = M(num_elem, blocksizes[i]);
+            }
+            return ForwardBlock(index, val, jac);
+        }
+
+        static ForwardBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
+        {
+
+            std::vector<M> jac;
+            const int num_elem = val.size();
+            const int num_blocks = blocksizes.size();
+            // First, set all jacobian blocks to zero...
+            jac.resize(num_blocks);
+            for (int i = 0; i < num_blocks; ++i) {
+                jac[i] = M(num_elem, blocksizes[i]);
+            }
+            // ... then set the one corrresponding to this variable to identity.
+            jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1));
+            for (typename M::Index row = 0; row < val.size(); ++row) {
+                jac[index].insert(row, row) = Scalar(1.0);
+            }
+            return ForwardBlock(index, val, jac);
+        }
+
+        static ForwardBlock function(const int index, const V& val, const std::vector<M>& jac)
+        {
+            return ForwardBlock(index, val, jac);
+        }
+
+#if 0
+        /// Operators.
+        ForwardBlock operator+(const ForwardBlock& rhs)
+        {
+            return function(val_ + rhs.val_, jac_ + rhs.jac_);
+        }
+
+        /// Operators.
+        ForwardBlock operator-(const ForwardBlock& rhs)
+        {
+            return function(val_ - rhs.val_, jac_ - rhs.jac_);
+        }
+
+        /// Operators.
+        ForwardBlock operator*(const ForwardBlock& rhs)
+        {
+            typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
+            D D1 = val_.matrix().asDiagonal();
+            D D2 = rhs.val_.matrix().asDiagonal();
+            return function(val_ * rhs.val_, D2*jac_ + D1*rhs.jac_);
+        }
+
+        /// Operators.
+        ForwardBlock operator/(const ForwardBlock& rhs)
+        {
+            typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
+            D D1 = val_.matrix().asDiagonal();
+            D D2 = rhs.val_.matrix().asDiagonal();
+            D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal();
+            return function(val_ / rhs.val_, D3 * (D2*jac_ - D1*rhs.jac_));
+        }
+#endif
+        /// I/O.
+        template <class Ostream>
+        Ostream&
+        print(Ostream& os) const
+        {
+            int num_blocks = jac_.size();
+            os << "Variable index " << index_ << ":\nValue =\n" << val_ << "\n\nJacobian =\n";
+            for (int i = 0; i < num_blocks; ++i) {
+                os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n";
+            }
+            return os;
+        }
+
+    private:
+        ForwardBlock(const int index,
+                     const V& val,
+                     const std::vector<M>& jac)
+            : index_(index), val_(val), jac_(jac)
+        {
+            const int num_elem = val_.size();
+            const int num_blocks = jac_.size();
+            assert(index_ < num_blocks);
+            assert(num_elem == jac_[index].cols());
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(num_elem == jac_[block].rows());
+            }
+        }
+
+        int index_;
+        V val_;
+        std::vector<M> jac_;
+    };
+
+
+    template <class Ostream, typename Scalar>
+    Ostream&
+    operator<<(Ostream& os, const ForwardBlock<Scalar>& fw)
+    {
+        return fw.print(os);
+    }
+
+
+
+} // namespace Autodiff
+
+
+
+#endif // OPM_AUTODIFFBLOCK_HEADER_INCLUDED

test_block.cpp

@@ -0,0 +1,60 @@
+/*
+  Copyright 2013 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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "AutoDiffBlock.hpp"
+#include <iostream>
+
+int main()
+{
+    typedef AutoDiff::ForwardBlock<double> ADV;
+    std::vector<int> blocksizes = { 3, 1, 2 };
+    int num_blocks = blocksizes.size();
+    ADV::V v1(3);
+    v1 << 0.2, 1.2, 13.4;
+    ADV::V v2(3);
+    v2 << 1.0, 2.2, 3.4;
+    enum { FirstVar = 0, SecondVar = 1, ThirdVar = 2 };
+    ADV a = ADV::constant(FirstVar, v1, blocksizes);
+    ADV x = ADV::variable(FirstVar, v2, blocksizes);
+    std::vector<ADV::M> jacs(num_blocks);
+    for (int i = 0; i < num_blocks; ++i) {
+        jacs[i] = ADV::M(blocksizes[FirstVar], blocksizes[i]);
+        jacs[i].insert(0,0) = -1.0;
+    }
+    ADV f = ADV::function(FirstVar, v2, jacs);
+    std::cout << a << x << f;
+    /*
+    ADV xpx = x + x;
+    std::cout << xpx;
+    ADV xpxpa = x + x + a;
+    std::cout << xpxpa;
+
+    std::cout << xpxpa - xpx;
+
+    ADV sqx = x * x;
+
+    std::cout << sqx;
+
+    ADV sqxdx = sqx / x;
+
+    std::cout << sqxdx;
+
+    // std::cout << a << "\n\n" << x << "\n\n" << f << std::endl;
+    */
+}