]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add new class TensorProductMatrix
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 10 Jul 2017 19:00:27 +0000 (21:00 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 14 Jul 2017 12:29:09 +0000 (14:29 +0200)
Tests for new TensorProductMatrix class

include/deal.II/lac/tensor_product_matrix.h [new file with mode: 0644]
tests/lac/tensor_product_matrix_01.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_01.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_02.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_02.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_03.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_03.wih_lapack=true.output [new file with mode: 0644]

diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h
new file mode 100644 (file)
index 0000000..24b920c
--- /dev/null
@@ -0,0 +1,412 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__tensor_product_matrix_h
+#define dealii__tensor_product_matrix_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <typename> class Vector;
+template <typename> class FullMatrix;
+
+/**
+ * This is a special matrix class defined as the tensor product (or Kronecker
+ * product) of 1D matrices of the type
+ * @f{align*}{
+ * L &= A \otimes M + M \otimes A
+ * @f}
+ * in 2D and
+ * @f{align*}{
+ * L &= A \otimes M \otimes M + M \otimes A \otimes M + M \otimes M \otimes A
+ * @f}
+ * in 3D. The typical application setting is a discretization of the Laplacian
+ * $L$ on a Cartesian (axis-aligned) geometry, where it can be exactly
+ * represented by the Kronecker or tensor product of a 1D mass matrix $M$ and
+ * a 1D Laplace matrix $A$ in each dimension. The dimension of the resulting
+ * class is the product of the one-dimensional matrices.
+ *
+ * This class implements two basic operations, namely the usual multiplication
+ * by a vector and the inverse. For both operations, fast tensorial techniques
+ * can be applied that implement the operator evaluation in
+ * $\text{size}(M)^{d+1}$ arithmetic operations, considerably less than
+ * $\text{size}(M)^{2d}$ for the naive forward transformation and
+ * $\text{size}(M)^{3d}$ for setting up the inverse of $L$.
+ *
+ * Interestingly, the exact inverse of the matrix $L$ can be found through
+ * tensor products due to an article by <a
+ * href="http://dl.acm.org/citation.cfm?id=2716130">R. E. Lynch, J. R. Rice,
+ * D. H. Thomas, Direct solution of partial difference equations by tensor
+ * product methods, Numerische Mathematik 6, 185-199</a> from 1964,
+ * @f{align*}{
+ * L^{-1} &= S \otimes S (\Lambda \otimes I + I \otimes \Lambda)^{-1}
+ * S^\mathrm T \otimes S^\mathrm T,
+ * @f}
+ * where $S$ is the matrix of eigenvectors to the generalized eigenvalue problem
+ * @f{align*}{
+ * A s  &= \lambda M s,
+ * @f}
+ * and $\Lambda$ is the diagonal matrix representing the generalized
+ * eigenvalues $\lambda$. Note that the vectors $s$ are such that they
+ * simultaneously diagonalize $A$ and $M$, $S^{\mathrm T} A S = \Lambda$ and
+ * $S^{\mathrm T} B S = I$. This method of matrix inversion is called fast
+ * diagonalization method.
+ *
+ * This class requires LAPACK support.
+ *
+ * Note that this class allows for two modes of usage. The first is a use case
+ * with run time constants for the matrix dimensions that is achieved by
+ * setting the optional template parameter for the size to -1. The second mode
+ * of usage that is faster allows to set the template parameter as a compile
+ * time constant, giving significantly faster code in particular for small
+ * sizes of the matrix.
+ *
+ * @note This class uses a temporary array for storing intermediate results
+ * that is a class member. A mutex is used to protect access to this array and
+ * ensure correct results. If several threads run parallel instances of this
+ * class, it is recommended that each threads holds its own matrix version.
+ *
+ * @tparam dim Dimension of the problem. Currently, 1D, 2D, and 3D codes are
+ * implemented.
+ *
+ * @tparam Number Type of the underlying array elements. Note that the
+ * underlying LAPACK implementation supports only float and double numbers, so
+ * only these two types are currently supported.
+ *
+ * @tparam size Compile-time array lengths. By default at -1, which means that
+ * the run-time info stored in the matrices passed to the reinit()
+ * function is used.
+ *
+ * @author Martin Kronbichler, 2017
+ */
+template <int dim, typename Number, int size = -1>
+class TensorProductMatrixSymmetricSum
+{
+public:
+  /**
+   * Constructor.
+   */
+  TensorProductMatrixSymmetricSum();
+
+  /**
+   * Constructor that is equivalent to the previous constructor and
+   * immediately calling reinit().
+   */
+  TensorProductMatrixSymmetricSum(const FullMatrix<Number> &mass_matrix,
+                                  const FullMatrix<Number> &derivative_matrix);
+
+  /**
+   * Initializes the matrix to the given mass matrix $M$ and derivative matrix
+   * $A$. Note that the current implementation requires $M$ to be symmetric
+   * and positive definite and $A$ to be symmetric and invertible but not
+   * necessarily positive defininte.
+   */
+  void reinit (const FullMatrix<Number> &mass_matrix,
+               const FullMatrix<Number> &derivative_matrix);
+
+  /**
+   * Returns the number of rows of this matrix, given by the dim-th power of
+   * the size of the 1D matrices passed to the constructor.
+   */
+  unsigned int m() const;
+
+  /**
+   * Returns the number of columns of this matrix, given by the dim-th power
+   * of the size of the 1D matrices passed to the constructor.
+   */
+  unsigned int n() const;
+
+  /**
+   * Implements a matrix-vector product with the underlying matrix as
+   * described in the main documentation of this class.
+   */
+  void vmult (Vector<Number> &dst,
+              const Vector<Number> &src) const;
+
+  /**
+   * Implements a matrix-vector product with the underlying matrix as
+   * described in the main documentation of this class. Same as the other
+   * vmult() function, but operating on plain pointers rather than a vector
+   * (no check of array bounds possible).
+   */
+  void vmult (Number *dst,
+              const Number *src) const;
+
+  /**
+   * Implements a matrix-vector product with the underlying matrix as
+   * described in the main documentation of this class.
+   */
+  void apply_inverse (Vector<Number> &dst,
+                      const Vector<Number> &src) const;
+
+  /**
+   * Implements a matrix-vector product with the underlying matrix as
+   * described in the main documentation of this class. Same as the other
+   * apply_inverse() function, but operating on plain pointers rather than a
+   * vector (no check of array bounds possible).
+   */
+  void apply_inverse (Number *dst,
+                      const Number *src) const;
+
+private:
+  /**
+   * A copy of the @p mass_matrix object passed to the reinit() method.
+   */
+  FullMatrix<Number> mass_matrix;
+
+  /**
+   * A copy of the @p derivative_matrix object passed to the reinit() method.
+   */
+  FullMatrix<Number> derivative_matrix;
+
+  /**
+   * A vector containing the generalized eigenvalues of A s = lambda B s.
+   */
+  AlignedVector<Number> eigenvalues;
+
+  /**
+   * The matrix containing the generalized eigenvectors.
+   */
+  Table<2,Number> eigenvectors;
+
+  /**
+   * An array for temporary data.
+   */
+  mutable AlignedVector<Number> tmp_array;
+
+  /**
+   * A mutex that guards access to the array @p tmp_array.
+   */
+  mutable Threads::Mutex mutex;
+};
+
+
+/*----------------------- Inline functions ----------------------------------*/
+
+#ifndef DOXYGEN
+
+
+template <int dim, typename Number, int size>
+inline
+TensorProductMatrixSymmetricSum<dim,Number,size>
+::TensorProductMatrixSymmetricSum()
+{}
+
+
+
+template <int dim, typename Number, int size>
+inline
+TensorProductMatrixSymmetricSum<dim,Number,size>
+::TensorProductMatrixSymmetricSum(const FullMatrix<Number> &mass_matrix,
+                                  const FullMatrix<Number> &derivative_matrix)
+{
+  reinit(mass_matrix, derivative_matrix);
+}
+
+
+
+template <int dim, typename Number, int size>
+inline
+void
+TensorProductMatrixSymmetricSum<dim,Number,size>
+::reinit(const FullMatrix<Number> &mass_matrix,
+         const FullMatrix<Number> &derivative_matrix)
+{
+  Assert(size == -1 ||
+         (size > 0 && static_cast<unsigned int>(size) == mass_matrix.m()),
+         ExcDimensionMismatch(size, mass_matrix.m()));
+  AssertDimension(mass_matrix.m(), mass_matrix.n());
+  AssertDimension(mass_matrix.m(), derivative_matrix.m());
+  AssertDimension(mass_matrix.m(), derivative_matrix.n());
+
+  this->mass_matrix = mass_matrix;
+  this->derivative_matrix = derivative_matrix;
+
+  std::vector<Vector<Number> > eigenvecs(mass_matrix.m());
+  LAPACKFullMatrix<Number> mass_copy(mass_matrix.m(), mass_matrix.n());
+  LAPACKFullMatrix<Number> deriv_copy(derivative_matrix.m(), derivative_matrix.n());
+  mass_copy = mass_matrix;
+  deriv_copy = derivative_matrix;
+
+  deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy, eigenvecs);
+  AssertDimension(eigenvecs.size(), mass_matrix.m());
+  eigenvectors.reinit(mass_matrix.m(), mass_matrix.m());
+  for (unsigned int i=0; i<mass_matrix.m(); ++i)
+    for (unsigned int j=0; j<mass_matrix.n(); ++j)
+      eigenvectors(i,j) = eigenvecs[j][i];
+
+  eigenvalues.resize(mass_matrix.m());
+  for (unsigned int i=0; i<mass_matrix.m(); ++i)
+    eigenvalues[i] = deriv_copy.eigenvalue(i).real();
+}
+
+
+
+template <int dim, typename Number, int size>
+inline
+unsigned int
+TensorProductMatrixSymmetricSum<dim,Number,size>::m() const
+{
+  return Utilities::fixed_power<dim>(mass_matrix.m());
+}
+
+
+
+template <int dim, typename Number, int size>
+inline
+unsigned int
+TensorProductMatrixSymmetricSum<dim,Number,size>::n() const
+{
+  return Utilities::fixed_power<dim>(mass_matrix.n());
+}
+
+
+
+template <int dim, typename Number, int size>
+inline
+void
+TensorProductMatrixSymmetricSum<dim,Number,size>
+::vmult(Vector<Number> &dst,
+        const Vector<Number> &src) const
+{
+  AssertDimension(dst.size(), Utilities::fixed_power<dim>(eigenvalues.size()));
+  AssertDimension(src.size(), Utilities::fixed_power<dim>(eigenvalues.size()));
+  vmult(dst.begin(), src.begin());
+}
+
+
+
+template <int dim, typename Number, int size>
+inline
+void
+TensorProductMatrixSymmetricSum<dim,Number,size>
+::apply_inverse(Vector<Number> &dst,
+                const Vector<Number> &src) const
+{
+  AssertDimension(dst.size(), Utilities::fixed_power<dim>(eigenvalues.size()));
+  AssertDimension(src.size(), Utilities::fixed_power<dim>(eigenvalues.size()));
+  apply_inverse(dst.begin(), src.begin());
+}
+
+
+
+template <int dim, typename Number, int size>
+inline
+void
+TensorProductMatrixSymmetricSum<dim,Number,size>
+::vmult(Number *dst,
+        const Number *src) const
+{
+  Threads::Mutex::ScopedLock lock(this->mutex);
+  const unsigned int n = Utilities::fixed_power<dim>(size > 0 ? size : eigenvalues.size());
+  tmp_array.resize_fast(n*2);
+  const int kernel_size = size > 0 ? size-1 : -1;
+  internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size+1,Number>
+  eval(AlignedVector<Number>(), AlignedVector<Number>(),
+       AlignedVector<Number>(), mass_matrix.m()-1, mass_matrix.m());
+  const Number *A = &derivative_matrix(0,0);
+  const Number *M = &mass_matrix(0,0);
+  Number *t = tmp_array.begin();
+  if (dim == 1)
+    eval.template apply<0, true, false>(A, src, dst);
+  else if (dim == 2)
+    {
+      eval.template apply<0, true, false>(M, src, t);
+      eval.template apply<1, true, false>(A, t, dst);
+      eval.template apply<0, true, false>(A, src, t);
+      eval.template apply<1, true, true> (M, t, dst);
+    }
+  else if (dim == 3)
+    {
+      eval.template apply<0, true, false>(M, src, t+n);
+      eval.template apply<1, true, false>(M, t+n, t);
+      eval.template apply<2, true, false>(A, t, dst);
+      eval.template apply<1, true, false>(A, t+n, t);
+      eval.template apply<0, true, false>(A, src, t+n);
+      eval.template apply<1, true, true> (M, t+n, t);
+      eval.template apply<2, true, true> (M, t, dst);
+    }
+  else
+    AssertThrow(false, ExcNotImplemented());
+}
+
+
+
+template <int dim, typename Number, int size>
+inline
+void
+TensorProductMatrixSymmetricSum<dim,Number,size>
+::apply_inverse(Number *dst,
+                const Number *src) const
+{
+  Threads::Mutex::ScopedLock lock(this->mutex);
+  const unsigned int n = size > 0 ? size : eigenvalues.size();
+  tmp_array.resize_fast(Utilities::fixed_power<dim>(n));
+  const int kernel_size = size > 0 ? size-1 : -1;
+  internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size+1,Number>
+  eval(AlignedVector<Number>(), AlignedVector<Number>(),
+       AlignedVector<Number>(), mass_matrix.m()-1, mass_matrix.m());
+  const Number *S = &eigenvectors(0,0);
+  Number *t = tmp_array.begin();
+
+  switch (dim)
+    {
+    case 1:
+      eval.template apply<0, true, false> (S, src, t);
+      for (unsigned int i=0; i<n; ++i)
+        t[i] /= eigenvalues[i];
+      eval.template apply<0, false, false> (S, t, dst);
+      break;
+
+    case 2:
+      eval.template apply<0, true, false> (S, src, t);
+      eval.template apply<1, true, false> (S, t, dst);
+      for (unsigned int i=0, c=0; i<n; ++i)
+        for (unsigned int j=0; j<n; ++j, ++c)
+          dst[c] /= (eigenvalues[i] + eigenvalues[j]);
+      eval.template apply<1, false, false> (S, dst, t);
+      eval.template apply<0, false, false> (S, t, dst);
+      break;
+
+    case 3:
+      eval.template apply<0, true, false> (S, src, t);
+      eval.template apply<1, true, false> (S, t, dst);
+      eval.template apply<2, true, false> (S, dst, t);
+      for (unsigned int i=0, c=0; i<n; ++i)
+        for (unsigned int j=0; j<n; ++j)
+          for (unsigned int k=0; k<n; ++k, ++c)
+            t[c] /= (eigenvalues[i] + eigenvalues[j] + eigenvalues[k]);
+      eval.template apply<2, false, false> (S, t, dst);
+      eval.template apply<1, false, false> (S, dst, t);
+      eval.template apply<0, false, false> (S, t, dst);
+      break;
+
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
+#endif
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/tests/lac/tensor_product_matrix_01.cc b/tests/lac/tensor_product_matrix_01.cc
new file mode 100644 (file)
index 0000000..c618001
--- /dev/null
@@ -0,0 +1,103 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test non-templated path of TensorProductMatrix
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+template <int dim>
+void do_test(const unsigned int size)
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  FullMatrix<double> mass(size, size);
+  FullMatrix<double> laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      mass(i,i) = 2./3.;
+      if (i > 0)
+        mass(i,i-1) = 1./6.;
+      if (i<size-1)
+        mass(i,i+1) = 1./6.;
+      laplace(i,i) = 2.;
+      if (i > 0)
+        laplace(i,i-1) = -1.;
+      if (i < size-1)
+        laplace(i,i+1) = -1.;
+    }
+  TensorProductMatrixSymmetricSum<dim,double> mat;
+  mat.reinit(mass, laplace);
+  Vector<double> v1(mat.m()), v2(mat.m()), v3(mat.m());
+  for (unsigned int i=0; i<v1.size(); ++i)
+    v1(i) = (2*i+1)%23;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  v3 -= v1;
+  deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
+
+  FullMatrix<double> full(v1.size(), v1.size());
+  for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+    for (unsigned int j=0; j<(dim>1?size:1); ++j)
+      for (unsigned int k=0; k<size; ++k, ++c)
+        for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+          for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+            for (unsigned int kk=0; kk<size; ++kk, ++cc)
+              if (dim == 1)
+                full(c,cc) = laplace(k,kk);
+              else if (dim==2)
+                full(c,cc) = laplace(k,kk)*mass(j,jj) + laplace(j,jj)*mass(k,kk);
+              else if (dim==3)
+                full(c,cc) = (laplace(k,kk)*mass(j,jj) + laplace(j,jj)*mass(k,kk))*mass(i,ii)
+                             + laplace(i,ii)*mass(j,jj)*mass(k,kk);
+  full.vmult(v3, v1);
+  v3 -= v2;
+  deallog << "Verifiction of vmult: " << v3.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(v3, v1);
+  mat.apply_inverse(v2, v1);
+  v3 -= v2;
+  deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1e-11);
+
+  do_test<1>(1);
+  do_test<1>(2);
+  do_test<1>(5);
+  do_test<2>(1);
+  do_test<2>(2);
+  do_test<2>(5);
+  do_test<2>(11);
+  do_test<3>(1);
+  do_test<3>(2);
+  do_test<3>(3);
+  do_test<3>(7);
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_01.with_lapack=true.output b/tests/lac/tensor_product_matrix_01.with_lapack=true.output
new file mode 100644 (file)
index 0000000..bde376c
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
diff --git a/tests/lac/tensor_product_matrix_02.cc b/tests/lac/tensor_product_matrix_02.cc
new file mode 100644 (file)
index 0000000..e5c5fb0
--- /dev/null
@@ -0,0 +1,103 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test templated path of TensorProductMatrix
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+template <int dim, int size>
+void do_test()
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  FullMatrix<double> mass(size, size);
+  FullMatrix<double> laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      mass(i,i) = 2./3.;
+      if (i > 0)
+        mass(i,i-1) = 1./6.;
+      if (i<size-1)
+        mass(i,i+1) = 1./6.;
+      laplace(i,i) = 2.;
+      if (i > 0)
+        laplace(i,i-1) = -1.;
+      if (i < size-1)
+        laplace(i,i+1) = -1.;
+    }
+  TensorProductMatrixSymmetricSum<dim,double,size> mat;
+  mat.reinit(mass, laplace);
+  Vector<double> v1(mat.m()), v2(mat.m()), v3(mat.m());
+  for (unsigned int i=0; i<v1.size(); ++i)
+    v1(i) = (2*i+1)%23;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  v3 -= v1;
+  deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
+
+  FullMatrix<double> full(v1.size(), v1.size());
+  for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+    for (unsigned int j=0; j<(dim>1?size:1); ++j)
+      for (unsigned int k=0; k<size; ++k, ++c)
+        for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+          for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+            for (unsigned int kk=0; kk<size; ++kk, ++cc)
+              if (dim == 1)
+                full(c,cc) = laplace(k,kk);
+              else if (dim==2)
+                full(c,cc) = laplace(k,kk)*mass(j,jj) + laplace(j,jj)*mass(k,kk);
+              else if (dim==3)
+                full(c,cc) = (laplace(k,kk)*mass(j,jj) + laplace(j,jj)*mass(k,kk))*mass(i,ii)
+                             + laplace(i,ii)*mass(j,jj)*mass(k,kk);
+  full.vmult(v3, v1);
+  v3 -= v2;
+  deallog << "Verifiction of vmult: " << v3.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(v3, v1);
+  mat.apply_inverse(v2, v1);
+  v3 -= v2;
+  deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1e-11);
+
+  do_test<1,1>();
+  do_test<1,2>();
+  do_test<1,5>();
+  do_test<2,1>();
+  do_test<2,2>();
+  do_test<2,5>();
+  do_test<2,11>();
+  do_test<3,1>();
+  do_test<3,2>();
+  do_test<3,3>();
+  do_test<3,7>();
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_02.with_lapack=true.output b/tests/lac/tensor_product_matrix_02.with_lapack=true.output
new file mode 100644 (file)
index 0000000..bde376c
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
diff --git a/tests/lac/tensor_product_matrix_03.cc b/tests/lac/tensor_product_matrix_03.cc
new file mode 100644 (file)
index 0000000..e06ed51
--- /dev/null
@@ -0,0 +1,103 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test templated path of TensorProductMatrix with float instead of double
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+template <int dim, int size>
+void do_test()
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  FullMatrix<float> mass(size, size);
+  FullMatrix<float> laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      mass(i,i) = 2./3.;
+      if (i > 0)
+        mass(i,i-1) = 1./6.;
+      if (i<size-1)
+        mass(i,i+1) = 1./6.;
+      laplace(i,i) = 2.;
+      if (i > 0)
+        laplace(i,i-1) = -1.;
+      if (i < size-1)
+        laplace(i,i+1) = -1.;
+    }
+  TensorProductMatrixSymmetricSum<dim,float,size> mat;
+  mat.reinit(mass, laplace);
+  Vector<float> v1(mat.m()), v2(mat.m()), v3(mat.m());
+  for (unsigned int i=0; i<v1.size(); ++i)
+    v1(i) = (2*i+1)%23;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  v3 -= v1;
+  deallog << "Verification of vmult and inverse: " << (double)v3.linfty_norm() << std::endl;
+
+  FullMatrix<float> full(v1.size(), v1.size());
+  for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+    for (unsigned int j=0; j<(dim>1?size:1); ++j)
+      for (unsigned int k=0; k<size; ++k, ++c)
+        for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+          for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+            for (unsigned int kk=0; kk<size; ++kk, ++cc)
+              if (dim == 1)
+                full(c,cc) = laplace(k,kk);
+              else if (dim==2)
+                full(c,cc) = laplace(k,kk)*mass(j,jj) + laplace(j,jj)*mass(k,kk);
+              else if (dim==3)
+                full(c,cc) = (laplace(k,kk)*mass(j,jj) + laplace(j,jj)*mass(k,kk))*mass(i,ii)
+                             + laplace(i,ii)*mass(j,jj)*mass(k,kk);
+  full.vmult(v3, v1);
+  v3 -= v2;
+  deallog << "Verifiction of vmult: " << (double)v3.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(v3, v1);
+  mat.apply_inverse(v2, v1);
+  v3 -= v2;
+  deallog << "Verification of inverse: " << (double)v3.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1e-3);
+
+  do_test<1,1>();
+  do_test<1,2>();
+  do_test<1,5>();
+  do_test<2,1>();
+  do_test<2,2>();
+  do_test<2,5>();
+  do_test<2,11>();
+  do_test<3,1>();
+  do_test<3,2>();
+  do_test<3,3>();
+  do_test<3,7>();
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_03.wih_lapack=true.output b/tests/lac/tensor_product_matrix_03.wih_lapack=true.output
new file mode 100644 (file)
index 0000000..bde376c
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.