From e2d8b77ec537a6410ca108f8ba3ec25235182701 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Tue, 30 Apr 2013 11:48:04 +0200
Subject: [PATCH] Initial vector AD type experimentation.

---
 AutoDiffVec.hpp | 97 +++++++++++++++++++++++++++++++++++++++++++++++++
 test_vec.cpp    | 36 ++++++++++++++++++
 2 files changed, 133 insertions(+)
 create mode 100644 AutoDiffVec.hpp
 create mode 100644 test_vec.cpp

diff --git a/AutoDiffVec.hpp b/AutoDiffVec.hpp
new file mode 100644
index 000000000..d555ae84a
--- /dev/null
+++ b/AutoDiffVec.hpp
@@ -0,0 +1,97 @@
+/*
+  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_AUTODIFFVEC_HEADER_INCLUDED
+#define OPM_AUTODIFFVEC_HEADER_INCLUDED
+
+#include "AutoDiff.hpp"
+#include <eigen3/Eigen/Eigen>
+#include <eigen3/Eigen/Sparse>
+#include <iostream>
+
+namespace AutoDiff
+{
+
+    template <typename Scalar>
+    class ForwardVec
+    {
+    public:
+        typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
+        typedef Eigen::SparseMatrix<Scalar> M;
+
+        static ForwardVec constant(const V& val)
+        {
+            return ForwardVec(val);
+        }
+
+        static ForwardVec variable(const V& val)
+        {
+            ForwardVec ret(val);
+            for (typename M::Index row = 0; row < val.size(); ++row) {
+                ret.jac_.insert(row, row) = 1.0;
+            }
+            ret.jac_.makeCompressed();
+            return ret;
+        }
+
+        static ForwardVec function(const V& val, const M& jac)
+        {
+            return ForwardVec(val, jac);
+        }
+
+        explicit ForwardVec(const V& val)
+            : val_(val), jac_(val.size(), val.size())
+        {
+        }
+
+        ForwardVec(const V& val, const M& jac)
+        : val_(val), jac_(jac)
+        {
+        }
+
+        template <class Ostream>
+        Ostream&
+        print(Ostream& os) const
+        {
+            os << "val =\n" << val_ << "\n\njac =\n" << jac_ << "\n";
+
+            return os;
+        }
+
+        
+    private:
+        V val_;
+        M jac_;
+    };
+
+
+    template <class Ostream, typename Scalar>
+    Ostream&
+    operator<<(Ostream& os, const ForwardVec<Scalar>& fw)
+    {
+        return fw.print(os);
+    }
+
+        
+
+} // namespace Autodiff
+
+
+
+#endif // OPM_AUTODIFFVEC_HEADER_INCLUDED
diff --git a/test_vec.cpp b/test_vec.cpp
new file mode 100644
index 000000000..6cd3f2287
--- /dev/null
+++ b/test_vec.cpp
@@ -0,0 +1,36 @@
+/*
+  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 "AutoDiffVec.hpp"
+#include <iostream>
+
+int main()
+{
+    typedef AutoDiff::ForwardVec<double> ADV;
+    ADV::V v(3);
+    v << 1.0, 2.2, 3.4;
+    std::cout << v << std::endl;
+    ADV a = ADV::constant(v);
+    ADV x = ADV::variable(v);
+    ADV::M jac(3,3);
+    jac.insert(0, 0) = 1.0;
+    ADV f = ADV::function(v, jac);
+
+    std::cout << a << "\n\n" << x << "\n\n" << f << std::endl;
+}