modified: defaults.cmake

modified:   ../examples/CMakeLists.txt
	new file:   ../fem/adnonlininteg.cpp
	new file:   ../fem/adnonlininteg.hpp
	modified:   ../linalg/CMakeLists.txt
	new file:   ../linalg/taddensemat.hpp
	new file:   ../linalg/tadvector.hpp

config/defaults.cmake

@@ -48,6 +48,8 @@ option(MFEM_USE_CUDA "Enable CUDA" OFF)
 option(MFEM_USE_OCCA "Enable OCCA" OFF)
 option(MFEM_USE_RAJA "Enable RAJA" OFF)
 option(MFEM_USE_CEED "Enable CEED" OFF)
+option(MFEM_USE_ADEPT "Enable AD using ADEPT"   OFF)
+option(MFEM_USE_CODIPACK "Enable AD using CoDiPack" OFF)
 
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
@@ -179,6 +181,13 @@ set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")
 set(LAPACK_LIBRARIES "" CACHE STRING "The LAPACK library.")
 
+set(ADEPT_INCLUDE_DIRS "${MFEM_DIR}/../adept-1.1op/include" CACHE STRING "Path to ADEPT headers.")
+set(ADEPT_LIBRARIES    "-L${MFEM_DIR}/../adept-1.1op/lib -ladept" CACHE STRING "The ADEPT library.")
+
+set(CODIPACK_INCLUDE_DIRS "${MFEM_DIR}/../CoDiPack/include" CACHE STRING "Path to CoDiPack headers.")
+
+
+
 # Some useful variables:
 set(CMAKE_SKIP_PREPROCESSED_SOURCE_RULES ON) # Skip *.i rules
 set(CMAKE_SKIP_ASSEMBLY_SOURCE_RULES  ON)    # Skip *.s rules

examples/CMakeLists.txt

@@ -121,3 +121,6 @@ if (MFEM_USE_HIOP)
   add_subdirectory(hiop)
 endif()
 
+if (MFEM_USE_ADEPT)
+  add_subdirectory(ad)
+endif()

fem/adnonlininteg.cpp

@@ -0,0 +1,33 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "fem.hpp"
+#include "../general/forall.hpp"
+#include "adnonlininteg.hpp"
+
+namespace mfem
+{
+
+ADNonlinearFormIntegrator::ADNonlinearFormIntegrator()
+{
+    
+}
+
+ADNonlinearFormIntegrator::~ADNonlinearFormIntegrator()
+{
+    
+}
+    
+    
+    
+    
+}
+

fem/adnonlininteg.hpp

@@ -0,0 +1,78 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+
+
+#ifndef MFEM_ADNONLININTEG
+#define MFEM_ADNONLININTEG
+
+#include "../config/config.hpp"
+#include "fe.hpp"
+#include "coefficient.hpp"
+#include "fespace.hpp"
+#include "nonlininteg.hpp"
+#include "tadvectro.hpp"
+#include "taddensmat.hpp"
+
+namespace mfem
+{
+
+#if define(MFEM_USE_ADEPT)||define(MFEM_USE_CODIPACK)
+/** The abstract base class ADNonlinearFormIntegrator is
+    a generalization of the NonlinearFormIntegrator class suitable
+    for algorithmic differentiation.   
+    All derived classes must implement ADAssembleElementVector(...); 
+    and ADGetElementEnergy(...); */
+class ADNonlinearFormIntegrator: public NonlinearFormIntegrator
+{
+protected:
+#ifdef  MFEM_USE_ADEPT
+     
+#elseif MFEM_USE_CODIPACK
+
+#endif    
+public:
+   ADNonlinearFormIntegrator();
+   virtual ~ADNonlinearFormIntegrator();
+ 
+   /// Methods called by the AD routines
+   virtual void ADGetElementEnergy(const mfem::FiniteElement & el,
+                                   mfem::ElementTransformation & Tr,
+                                   const mfem::TADVector<adouble> & elfun);
+   
+   virtual void ADAssembleElementVector(const mfem::FiniteElement & el,
+                                    mfem::ElementTransformation & Tr,
+                                    const mfem::TADVector<adouble> & elfun, 
+                                    mfem::TADVector<adouble> &elvec);
+   
+   
+   
+   /// Perform the local action of the NonlinearFormIntegrator
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &Tr,
+                                      const Vector &elfun, Vector &elvect) override;
+   
+   
+   virtual void AssembleElementGrad(const mfem::FiniteElement & el,
+                                    mfem::ElementTransformation & Tr,
+                                    const mfem::Vector & elfun, 
+                                    mfem::DenseMatrix & elmat) override;
+   
+   virtual double GetElementEnergy(const mfem::FiniteElement & el,
+                                   mfem::ElementTransformation & Tr,
+                                   const mfem::Vector & elfun) override;
+   
+};
+    
+    
+#endif    
+
+}

linalg/CMakeLists.txt

@@ -45,6 +45,8 @@ list(APPEND HDRS
   tmatrix.hpp
   ttensor.hpp
   vector.hpp
+  tadvector.hpp
+  taddensemat.hpp 
   )
 
 if (MFEM_USE_MPI)

linalg/taddensemat.hpp

@@ -0,0 +1,532 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+
+#ifndef TADDENSEMATRIX_H
+#define TADDENSEMATRIX_H
+
+
+#include "../config/config.hpp"
+#include "../general/globals.hpp"
+
+#include "tadvector.hpp"
+#include "densemat.hpp"
+
+
+
+namespace mfem
+{
+
+template<typename dtype>
+class TADDenseMatrix
+{
+private:
+   int height; ///< Dimension of the output / number of rows in the matrix.
+   int width;  ///< Dimension of the input / number of columns in the matrix.
+   dtype *data;
+   int capacity; // zero or negative capacity means we do not own the data.
+
+
+public:
+   /// Get the height (size of output) of the Operator. Synonym with NumRows().
+   inline int Height() const { return height; }
+   /** @brief Get the number of rows (size of output) of the Operator. Synonym
+       with Height(). */
+   inline int NumRows() const { return height; }
+
+   /// Get the width (size of input) of the Operator. Synonym with NumCols().
+   inline int Width() const { return width; }
+   /** @brief Get the number of columns (size of input) of the Operator. Synonym
+       with Width(). */
+   inline int NumCols() const { return width; }
+
+   /** Default constructor for TADDenseMatrix.
+      Sets data = NULL and height = width = 0. */
+   TADDenseMatrix()
+   {
+      data=nullptr;
+      capacity=0;
+      height=0;
+      width=0;
+   }
+
+   /// Copy constructor
+   template<typename idtype>
+   TADDenseMatrix(const TADDenseMatrix<idtype> &m)
+   {
+      height=m.GetHeight();
+      width=m.GetWidth();
+      const int hw = height * width;
+      if (hw > 0)
+      {
+         idtype* mdata=m.Data();
+         MFEM_ASSERT(mdata, "invalid source matrix");
+         data = new dtype[hw];
+         capacity = hw;
+         for (int i=0; i<hw; i++)
+         {
+            data[i]=mdata[i];
+         }
+      }
+      else
+      {
+         data = nullptr;
+         capacity = 0;
+         width=0;
+         height=0;
+      }
+   }
+
+   TADDenseMatrix(const DenseMatrix &m)
+   {
+      height=m.GetHeight();
+      width=m.GetWidth();
+      const int hw = height * width;
+      if (hw > 0)
+      {
+         double* mdata=m.Data();
+         MFEM_ASSERT(mdata, "invalid source matrix");
+         data = new dtype[hw];
+         capacity = hw;
+         for (int i=0; i<hw; i++)
+         {
+            data[i]=mdata[i];
+         }
+      }
+      else
+      {
+         data = nullptr;
+         capacity = 0;
+         width=0;
+         height=0;
+      }
+   }
+
+
+   /// Creates square matrix of size s.
+   explicit TADDenseMatrix(int s)
+   {
+      MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
+      height=s;
+      width=s;
+      capacity = s*s;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity](); // init with zeroes
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+   /// Creates rectangular matrix of size m x n.
+   TADDenseMatrix(int m, int n)
+   {
+      MFEM_ASSERT(m >= 0 && n >= 0,
+                  "invalid DenseMatrix size: " << m << " x " << n);
+      height=m;
+      width=n;
+      capacity = m*n;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity](); // init with zeroes
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+   TADDenseMatrix(const TADDenseMatrix<dtype> &mat, char ch)
+   {
+      height=mat.Width();
+      width=mat.Height();
+      capacity = height*width;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity];
+
+         for (int i = 0; i < height; i++)
+         {
+            for (int j = 0; j < width; j++)
+            {
+               (*this)(i,j) = mat(j,i);
+            }
+         }
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+
+   /// Change the size of the DenseMatrix to s x s.
+   void SetSize(int s) { SetSize(s, s); }
+
+   /// Change the size of the DenseMatrix to h x w.
+   void SetSize(int h, int w)
+   {
+      MFEM_ASSERT(h >= 0 && w >= 0,
+                  "invalid DenseMatrix size: " << h << " x " << w);
+      if (Height() == h && Width() == w)
+      {
+         return;
+      }
+      height = h;
+      width = w;
+      const int hw = h*w;
+      if (hw > std::abs(capacity))
+      {
+         if (capacity > 0)
+         {
+            delete [] data;
+         }
+         capacity = hw;
+         data = new dtype[hw](); // init with zeroes
+      }
+   }
+
+   /// Returns the matrix data array.
+   inline dtype *Data() const { return data; }
+   /// Returns the matrix data array.
+   inline dtype *GetData() const { return data; }
+
+   inline bool OwnsData() const { return (capacity > 0); }
+
+   /// Returns reference to a_{ij}.
+   dtype& operator()(int i, int j)
+   {
+      MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+      return data[i+j*height];
+   }
+
+   const dtype& operator()(int i, int j) const
+   {
+      MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+      return data[i+j*height];
+   }
+
+
+   dtype& Elem(int i, int j)
+   {
+      return (*this)(i,j);
+   }
+
+   const dtype& Elem(int i, int j) const
+   {
+      return (*this)(i,j);
+   }
+
+   void Mult(const dtype *x, dtype *y) const
+   {
+      if (width == 0)
+      {
+         for (int row = 0; row < height; row++)
+         {
+            y[row] = 0.0;
+         }
+         return;
+      }
+      dtype *d_col = data;
+      dtype x_col = x[0];
+      for (int row = 0; row < height; row++)
+      {
+         y[row] = x_col*d_col[row];
+      }
+      d_col += height;
+      for (int col = 1; col < width; col++)
+      {
+         x_col = x[col];
+         for (int row = 0; row < height; row++)
+         {
+            y[row] += x_col*d_col[row];
+         }
+         d_col += height;
+      }
+   }
+
+   void Mult(const TADVector<dtype> &x, TADVector<dtype> &y) const
+   {
+      MFEM_ASSERT(height == y.Size() && width == x.Size(),
+                  "incompatible dimensions");
+
+      Mult((const dtype *)x, (dtype *)y);
+   }
+
+   dtype operator *(const TADDenseMatrix<dtype> &m) const
+   {
+      MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
+                  "incompatible dimensions");
+
+      const int hw = height * width;
+      dtype a = 0.0;
+      for (int i = 0; i < hw; i++)
+      {
+         a += data[i] * m.data[i];
+      }
+
+      return a;
+   }
+
+   void MultTranspose(const dtype *x, dtype *y) const
+   {
+      dtype *d_col = data;
+      for (int col = 0; col < width; col++)
+      {
+         double y_col = 0.0;
+         for (int row = 0; row < height; row++)
+         {
+            y_col += x[row]*d_col[row];
+         }
+         y[col] = y_col;
+         d_col += height;
+      }
+   }
+
+   void MultTranspose(const  TADVector<dtype> &x,  TADVector<dtype> &y) const
+   {
+      MFEM_ASSERT(height == x.Size() && width == y.Size(),
+                  "incompatible dimensions");
+
+      MultTranspose((const dtype *)x, (dtype *)y);
+   }
+
+
+   void Randomize(int seed)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double)(RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int)time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned)seed);
+
+      for (int i = 0; i < capacity; i++)
+      {
+         data[i] = (dtype)(std::abs(rand()/max));
+      }
+   }
+
+   void RandomizeDiag(int seed)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double)(RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int)time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned)seed);
+
+      for (int i = 0; i < std::min(height,width); i++)
+      {
+         Elem(i,i) = (dtype)(std::abs(rand()/max));
+      }
+   }
+
+
+   /// Creates n x n diagonal matrix with diagonal elements c
+   void Diag(dtype c, int n)
+   {
+      SetSize(n);
+
+      const int N = n*n;
+      for (int i = 0; i < N; i++)
+      {
+         data[i] = (dtype)0.0;
+      }
+      for (int i = 0; i < n; i++)
+      {
+         data[i*(n+1)] = c;
+      }
+   }
+   /// Creates n x n diagonal matrix with diagonal given by diag
+   template<typename itype>
+   void Diag(itype *diag, int n)
+   {
+      SetSize(n);
+
+      int i, N = n*n;
+      for (i = 0; i < N; i++)
+      {
+         data[i] = 0.0;
+      }
+      for (i = 0; i < n; i++)
+      {
+         data[i*(n+1)] = (dtype) diag[i];
+      }
+   }
+
+   /// (*this) = (*this)^t
+   void Transpose()
+   {
+      int i, j;
+      dtype t;
+
+      if (Width() == Height())
+      {
+         for (i = 0; i < Height(); i++)
+            for (j = i+1; j < Width(); j++)
+            {
+               t = (*this)(i,j);
+               (*this)(i,j) = (*this)(j,i);
+               (*this)(j,i) = t;
+            }
+      }
+      else
+      {
+         TADDenseMatrix<dtype> T(*this,'t');
+         (*this) = T;
+      }
+   }
+   /// (*this) = A^t
+   template<typename itype>
+   void Transpose(const TADDenseMatrix<itype> &A)
+   {
+      SetSize(A.Width(),A.Height());
+
+      for (int i = 0; i < Height(); i++)
+         for (int j = 0; j < Width(); j++)
+         {
+            (*this)(i,j) = (dtype) A(j,i);
+         }
+   }
+
+
+   /// (*this) = 1/2 ((*this) + (*this)^t)
+   void Symmetrize()
+   {
+#ifdef MFEM_DEBUG
+      if (Width() != Height())
+      {
+         mfem_error("DenseMatrix::Symmetrize() : not a square matrix!");
+      }
+#endif
+
+      for (int i = 0; i < Height(); i++)
+         for (int j = 0; j < i; j++)
+         {
+            dtype a = 0.5 * ((*this)(i,j) + (*this)(j,i));
+            (*this)(j,i) = (*this)(i,j) = a;
+         }
+   }
+
+   void Lump()
+   {
+      for (int i = 0; i < Height(); i++)
+      {
+         dtype L = 0.0;
+         for (int j = 0; j < Width(); j++)
+         {
+            L += (*this)(i, j);
+            (*this)(i, j) = (dtype) 0.0;
+         }
+         (*this)(i, i) = L;
+      }
+   }
+
+
+
+
+};
+
+
+
+template<typename dtype>
+void CalcAdjugate(const TADDenseMatrix<dtype> &a, TADDenseMatrix<dtype> &adja)
+{
+#ifdef MFEM_DEBUG
+   if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
+   {
+      mfem_error("CalcAdjugate(...)");
+   }
+   if (a.Width() != adja.Height() || a.Height() != adja.Width())
+   {
+      mfem_error("CalcAdjugate(...)");
+   }
+#endif
+
+   if (a.Width() < a.Height())
+   {
+      const dtype *d = a.Data();
+      dtype *ad = adja.Data();
+      if (a.Width() == 1)
+      {
+         // N x 1, N = 2,3
+         ad[0] = d[0];
+         ad[1] = d[1];
+         if (a.Height() == 3)
+         {
+            ad[2] = d[2];
+         }
+      }
+      else
+      {
+         // 3 x 2
+         double e, g, f;
+         e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+         g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
+         f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
+
+         ad[0] = d[0]*g - d[3]*f;
+         ad[1] = d[3]*e - d[0]*f;
+         ad[2] = d[1]*g - d[4]*f;
+         ad[3] = d[4]*e - d[1]*f;
+         ad[4] = d[2]*g - d[5]*f;
+         ad[5] = d[5]*e - d[2]*f;
+      }
+      return;
+   }
+
+   if (a.Width() == 1)
+   {
+      adja(0,0) = (dtype)1.0;
+   }
+   else if (a.Width() == 2)
+   {
+      adja(0,0) =  a(1,1);
+      adja(0,1) = -a(0,1);
+      adja(1,0) = -a(1,0);
+      adja(1,1) =  a(0,0);
+   }
+   else
+   {
+      adja(0,0) = a(1,1)*a(2,2)-a(1,2)*a(2,1);
+      adja(0,1) = a(0,2)*a(2,1)-a(0,1)*a(2,2);
+      adja(0,2) = a(0,1)*a(1,2)-a(0,2)*a(1,1);
+
+      adja(1,0) = a(1,2)*a(2,0)-a(1,0)*a(2,2);
+      adja(1,1) = a(0,0)*a(2,2)-a(0,2)*a(2,0);
+      adja(1,2) = a(0,2)*a(1,0)-a(0,0)*a(1,2);
+
+      adja(2,0) = a(1,0)*a(2,1)-a(1,1)*a(2,0);
+      adja(2,1) = a(0,1)*a(2,0)-a(0,0)*a(2,1);
+      adja(2,2) = a(0,0)*a(1,1)-a(0,1)*a(1,0);
+   }
+}
+
+
+
+
+
+
+
+}
+
+#endif

linalg/tadvector.hpp

@@ -0,0 +1,687 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_TADVECTOR
+#define MFEM_TADVECTOR
+
+#include "../general/mem_manager.hpp"
+#include "vector.hpp"
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+#if defined(_MSC_VER) && (_MSC_VER < 1800)
+#include <float.h>
+#define isfinite _finite
+#endif
+
+
+namespace mfem
+{
+
+
+/// Vector data type.
+template<typename dtype>
+class TADVector
+{
+protected:
+
+   Memory<dtype> data;
+   int size;
+
+public:
+
+   /// Default constructor for Vector. Sets size = 0 and data = NULL.
+   TADVector() { data.Reset(); size = 0; }
+
+   /// Copy constructor. Allocates a new data array and copies the data.
+   TADVector(const TADVector<dtype> &v)
+   {
+      const int s = v.Size();
+      if (s > 0)
+      {
+         size = s;
+         data.New(s);
+         for (int i=0; i<s; i++)
+         {
+            data[i]=v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         data.Reset();
+      }
+   }
+
+   TADVector(const Vector &v)
+   {
+      const int s = v.Size();
+      if (s > 0)
+      {
+         size = s;
+         data.New(s);
+         for (int i=0; i<s; i++)
+         {
+            data[i]=v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         data.Reset();
+      }
+   }
+
+   /// @brief Creates vector of size s.
+   /// @warning Entries are not initialized to zero!
+   explicit TADVector(int s)
+   {
+      if (s > 0)
+      {
+         size = s;
+         data.New(s);
+      }
+      else
+      {
+         size = 0;
+         data.Reset();
+      }
+   }
+
+   /// Creates a vector referencing an array of doubles, owned by someone else.
+   /** The pointer @a _data can be NULL. The data array can be replaced later
+       with SetData(). */
+   TADVector(dtype *_data, int _size)
+   { data.Wrap(_data, _size, false); size = _size; }
+
+   /// Create a Vector of size @a size_ using MemoryType @a mt.
+   TADVector(int size_, MemoryType mt)
+      : data(size_, mt), size(size_) { }
+
+   /// Enable execution of Vector operations using the mfem::Device.
+   /** The default is to use Backend::CPU (serial execution on each MPI rank),
+       regardless of the mfem::Device configuration.
+
+       When appropriate, MFEM functions and class methods will enable the use
+       of the mfem::Device for their Vector parameters.
+
+       Some derived classes, e.g. GridFunction, enable the use of the
+       mfem::Device by default. */
+   void UseDevice(bool use_dev) const { data.UseDevice(use_dev); }
+
+   /// Return the device flag of the Memory object used by the Vector
+   bool UseDevice() const { return data.UseDevice(); }
+
+   /// Reads a vector from multiple files
+   void Load(std::istream ** in, int np, int * dim)
+   {
+      int i, j, s;
+
+      s = 0;
+      for (i = 0; i < np; i++)
+      {
+         s += dim[i];
+      }
+      SetSize(s);
+
+      int p = 0;
+      double tmpd;
+      for (i = 0; i < np; i++)
+      {
+         for (j = 0; j < dim[i]; j++)
+         {
+            *in[i] >> tmpd;
+            data[p++]=dtype(tmpd);
+         }
+      }
+   }
+
+   /// Load a vector from an input stream.
+   void Load(std::istream &in, int Size)
+   {
+      SetSize(Size);
+      double tmpd;
+      for (int i = 0; i < size; i++)
+      {
+         in >> tmpd;
+         data[i]=dtype(tmpd);
+      }
+   }
+
+   /// Load a vector from an input stream, reading the size from the stream.
+   void Load(std::istream &in) { int s; in >> s; Load(in, s); }
+
+   /// @brief Resize the vector to size @a s.
+   /** If the new size is less than or equal to Capacity() then the internal
+       data array remains the same. Otherwise, the old array is deleted, if
+       owned, and a new array of size @a s is allocated without copying the
+       previous content of the Vector.
+       @warning In the second case above (new size greater than current one),
+       the vector will allocate new data array, even if it did not own the
+       original data! Also, new entries are not initialized! */
+   void SetSize(int s)
+   {
+      if (s == size)
+      {
+         return;
+      }
+      if (s <= data.Capacity())
+      {
+         size = s;
+         return;
+      }
+      // preserve a valid MemoryType and device flag
+      const MemoryType mt = data.GetMemoryType();
+      const bool use_dev = data.UseDevice();
+      data.Delete();
+      size = s;
+      data.New(s, mt);
+      data.UseDevice(use_dev);
+   }
+
+   /// Resize the vector to size @a s using MemoryType @a mt.
+   void SetSize(int s, MemoryType mt)
+   {
+      if (mt == data.GetMemoryType())
+      {
+         if (s == size)
+         {
+            return;
+         }
+         if (s <= data.Capacity())
+         {
+            size = s;
+            return;
+         }
+      }
+      const bool use_dev = data.UseDevice();
+      data.Delete();
+      if (s > 0)
+      {
+         data.New(s, mt);
+         size = s;
+      }
+      else
+      {
+         data.Reset();
+         size = 0;
+      }
+      data.UseDevice(use_dev);
+   }
+
+   /// Set the Vector data.
+   /// @warning This method should be called only when OwnsData() is false.
+   void SetData(dtype *d) { data.Wrap(d, data.Capacity(), false); }
+
+   /// Set the Vector data and size.
+   /** The Vector does not assume ownership of the new data. The new size is
+       also used as the new Capacity().
+       @warning This method should be called only when OwnsData() is false.
+       @sa NewDataAndSize(). */
+   void SetDataAndSize(dtype *d, int s)
+   { data.Wrap(d, s, false); size = s; }
+
+   /// Set the Vector data and size, deleting the old data, if owned.
+   /** The Vector does not assume ownership of the new data. The new size is
+       also used as the new Capacity().
+       @sa SetDataAndSize(). */
+   void NewDataAndSize(dtype *d, int s)
+   {
+      data.Delete();
+      SetDataAndSize(d, s);
+   }
+
+   /// Reset the Vector to use the given external Memory @a mem and size @a s.
+   /** If @a own_mem is false, the Vector will not own any of the pointers of
+       @a mem.
+       @sa NewDataAndSize(). */
+   void NewMemoryAndSize(const Memory<dtype> &mem, int s, bool own_mem)
+   {
+      data.Delete();
+      size = s;
+      data = mem;
+      if (!own_mem) { data.ClearOwnerFlags(); }
+
+   }
+
+   /// Reset the Vector to be a reference to a sub-vector of @a base.
+   inline void MakeRef(TADVector<dtype> &base, int offset, int size_)
+   {
+      data.Delete();
+      size = size_;
+      data.MakeAlias(base.GetMemory(), offset, size_);
+   }
+
+   /** @brief Reset the Vector to be a reference to a sub-vector of @a base
+       without changing its current size. */
+   inline void MakeRef(TADVector<dtype> &base, int offset)
+   {
+      data.Delete();
+      data.MakeAlias(base.GetMemory(), offset, size);
+   }
+
+   /// Set the Vector data (host pointer) ownership flag.
+   inline void MakeDataOwner() const { data.SetHostPtrOwner(true); }
+
+   /// Destroy a vector
+   void Destroy()
+   {
+      data.Delete();
+      size = 0;
+      data.Reset();
+   }
+
+   /// Returns the size of the vector.
+   inline int Size() const { return size; }
+
+   /// Return the size of the currently allocated data array.
+   /** It is always true that Capacity() >= Size(). */
+   inline int Capacity() const { return data.Capacity(); }
+
+   /// Return a pointer to the beginning of the Vector data.
+   /** @warning This method should be used with caution as it gives write access
+       to the data of const-qualified Vector%s. */
+   inline dtype *GetData() const
+   { return const_cast<dtype*>((const dtype*)data); }
+
+   /// Conversion to `double *`.
+   /** @note This conversion function makes it possible to use [] for indexing
+       in addition to the overloaded operator()(int). */
+   inline operator dtype *() { return data; }
+
+   /// Conversion to `const double *`.
+   /** @note This conversion function makes it possible to use [] for indexing
+       in addition to the overloaded operator()(int). */
+   inline operator const dtype *() const { return data; }
+
+   /// Return a reference to the Memory object used by the Vector.
+   Memory<dtype> &GetMemory() { return data; }
+
+   /** @brief Return a reference to the Memory object used by the Vector, const
+       version. */
+   const Memory<dtype> &GetMemory() const { return data; }
+
+   /// Update the memory location of the vector to match @a v.
+   void SyncMemory(const TADVector<dtype> &v) { GetMemory().Sync(v.GetMemory()); }
+
+   /// Update the alias memory location of the vector to match @a v.
+   void SyncAliasMemory(const TADVector<dtype> &v)
+   { GetMemory().SyncAlias(v.GetMemory(),Size()); }
+
+   /// Read the Vector data (host pointer) ownership flag.
+   inline bool OwnsData() const { return data.OwnsHostPtr(); }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline void StealData(dtype **p)
+   { *p = data; data.Reset(); size = 0; }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline dtype *StealData() { dtype *p; StealData(&p); return p; }
+
+   /// Access Vector entries. Index i = 0 .. size-1.
+   dtype &Elem(int i)
+   {
+      return operator()(i);
+   }
+   /// Read only access to Vector entries. Index i = 0 .. size-1.
+   const double &Elem(int i) const
+   {
+      return operator()(i);
+   }
+
+   /// Access Vector entries using () for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline double &operator()(int i)
+   {
+      MFEM_ASSERT(data && i >= 0 && i < size,
+                  "index [" << i << "] is out of range [0," << size << ")");
+
+      return data[i];
+   }
+
+   /// Read only access to Vector entries using () for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline const double &operator()(int i) const
+   {
+      MFEM_ASSERT(data && i >= 0 && i < size,
+                  "index [" << i << "] is out of range [0," << size << ")");
+
+      return data[i];
+   }
+
+   /// Dot product with a `dtype *` array.
+   dtype operator*(const dtype *v) const
+   {
+      dtype dot = 0.0;
+#ifdef MFEM_USE_LEGACY_OPENMP
+      #pragma omp parallel for reduction(+:dot)
+#endif
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+
+   /// Return the inner-product.
+   dtype operator*(const TADVector<dtype> &v) const
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      dtype dot = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+   dtype operator*(const Vector &v) const
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      dtype dot = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+   /// Copy Size() entries from @a v.
+   TADVector<dtype> &operator=(const dtype *v)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=v[i];
+      }
+      return *this;
+   }
+
+   /// Copy assignment.
+   /** @note Defining this method overwrites the implicitly defined copy
+       assignemnt operator. */
+   TADVector<dtype> &operator=(const TADVector<dtype> &v)
+   {
+      SetSize(v.Size());
+      for (int i=0; i<size; i++)
+      {
+         data[i]=v[i];
+      }
+      return *this;
+   }
+
+   TADVector<dtype> &operator=(const Vector &v)
+   {
+      SetSize(v.Size());
+      for (int i=0; i<size; i++)
+      {
+         data[i]=v[i];
+      }
+      return *this;
+   }
+
+   /// Redefine '=' for vector = constant.
+   template<typename ivtype>
+   TADVector &operator=(ivtype value)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=value;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator*=(ivtype c)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]*c;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator/=(ivtype c)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]/c;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator-=(ivtype c)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]-c;
+      }
+      return *this;
+   }
+
+   TADVector &operator-=(const TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]-v[i];
+      }
+      return *this;
+   }
+
+   TADVector &operator+=(const TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]+v[i];
+      }
+      return *this;
+   }
+
+   /// (*this) += a * Va
+   template<typename ivtype, typename vtype>
+   TADVector &Add(const ivtype a, const vtype &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]+a*v[i];
+      }
+      return *this;
+   }
+
+   /// (*this) = a * x
+   template<typename ivtype, typename vtype>
+   TADVector &Set(const ivtype a, const vtype &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=a*v[i];
+      }
+      return *this;
+   }
+
+   template<typename vtype>
+   void SetVector(const vtype &v, int offset)
+   {
+      MFEM_ASSERT(v.Size() + offset <= size, "invalid sub-vector");
+      for (int i = 0; i < size; i++)
+      {
+         data[i+offset] = v[i];
+      }
+   }
+
+   /// (*this) = -(*this)
+   void Neg()
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i]=-data[i];
+      }
+   }
+
+   /// Swap the contents of two Vectors
+   inline void Swap(TADVector &other)
+   {
+      Swap(data, other.data);
+      Swap(size, other.size);
+   }
+
+   /// Set v = v1 + v2.
+   template<typename vtype1, typename vtype2>
+   friend void add(const vtype1 &v1, const vtype2 &v2, TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(v1.Size() == v.Size(), "incompatible Vectors!");
+      MFEM_ASSERT(v2.Size() == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<v.Size(); i++)
+      {
+         v[i]=v1[i]+v2[i];
+      }
+   }
+
+   /// Set v = v1 + alpha * v2.
+   template<typename vtype1, typename ivtype, typename vtype2>
+   friend void add(const vtype1 &v1, ivtype alpha, const vtype2 &v2,
+                   TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(v1.Size() == v.Size(), "incompatible Vectors!");
+      MFEM_ASSERT(v2.Size() == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<v.Size(); i++)
+      {
+         v[i]=v1[i]+alpha*v2[i];
+      }
+   }
+
+
+   /// Destroys vector.
+   ~TADVector()
+   {
+      data.Delete();
+   }
+
+
+   /// Prints vector to stream out.
+   void Print(std::ostream &out = mfem::out, int width = 8) const
+   {
+      if (!size) { return; }
+      data.Read(MemoryClass::HOST, size);
+      for (int i = 0; 1; )
+      {
+         out << data[i];
+         i++;
+         if (i == size)
+         {
+            break;
+         }
+         if ( i % width == 0 )
+         {
+            out << '\n';
+         }
+         else
+         {
+            out << ' ';
+         }
+      }
+      out << '\n';
+   }
+
+   /// Set random values in the vector.
+   void Randomize(int seed = 0)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double)(RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int)time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned)seed);
+
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = std::abs(rand()/max);
+      }
+   }
+   /// Returns the l2 norm of the vector.
+   dtype Norml2() const
+   {
+      // Scale entries of Vector on the fly, using algorithms from
+      // std::hypot() and LAPACK's drm2. This scaling ensures that the
+      // argument of each call to std::pow is <= 1 to avoid overflow.
+      if (0 == size)
+      {
+         return 0.0;
+      } // end if 0 == size
+
+      if (1 == size)
+      {
+         return std::abs(data[0]);
+      } // end if 1 == size
+
+      dtype scale = 0.0;
+      dtype sum = 0.0;
+
+      for (int i = 0; i < size; i++)
+      {
+         if (data[i] != 0.0)
+         {
+            const dtype absdata = abs(data[i]);
+            if (scale <= absdata)
+            {
+               const dtype sqr_arg = scale / absdata;
+               sum = 1.0 + sum * (sqr_arg * sqr_arg);
+               scale = absdata;
+               continue;
+            } // end if scale <= absdata
+            const dtype sqr_arg = absdata / scale;
+            sum += (sqr_arg * sqr_arg); // else scale > absdata
+         } // end if data[i] != 0
+      }
+      return scale * sqrt(sum);
+   }
+
+   /// Returns the l_infinity norm of the vector.
+   dtype Normlinf() const
+   {
+      dtype max = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         max = max(abs(data[i]), max);
+      }
+      return max;
+   }
+   /// Returns the l_1 norm of the vector.
+   dtype Norml1() const
+   {
+      dtype sum = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         sum += abs(data[i]);
+      }
+      return sum;
+   }
+
+
+};
+
+
+
+} // namespace mfem
+
+#endif
+