]> https://gitweb.dealii.org/ - dealii.git/commitdiff
removing file depending on upper and lower triangular
authorESeNonFossiIo <esenonfossiio@gmail.com>
Tue, 16 Jun 2015 19:57:23 +0000 (21:57 +0200)
committerESeNonFossiIo <esenonfossiio@gmail.com>
Thu, 25 Jun 2015 13:48:09 +0000 (15:48 +0200)
include/deal.II/lac/block_linear_operator.h
tests/lac/block_linear_operator_03.cc [new file with mode: 0644]
tests/lac/block_linear_operator_03.with_cxx11=on.output [moved from tests/lac/block_linear_operator_06.with_cxx11=on.output with 100% similarity]
tests/lac/block_linear_operator_04.cc [moved from tests/lac/block_linear_operator_09.cc with 96% similarity]
tests/lac/block_linear_operator_04.with_cxx11=on.output [moved from tests/lac/block_linear_operator_09.with_cxx11=on.output with 100% similarity]
tests/lac/block_linear_operator_05.cc [new file with mode: 0644]
tests/lac/block_linear_operator_05.with_cxx11=on.output [new file with mode: 0644]
tests/lac/block_linear_operator_06.cc [deleted file]

index 1a480a7d452919c1f3e670113b75cbc6acd2150a..64ff277b8f24fae65dc5306750ba6deb9fad387a 100644 (file)
@@ -322,352 +322,13 @@ block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename
 }
 
 
-/**
- * @relates LinearOperator
- *
- * A function that takes an array of arrays of LinearOperators @p block_matrix
- * and returns its associated lower triangular matrix operator
- * (diagonal is not included).
- *
- * @code
- *  a00 | a01 | a02                 |     |
- *  ---------------             ---------------
- *  a10 | a11 | a12     ->      a10 |     |
- *  ---------------             ---------------
- *  a20 | a21 | a22             a20 | a21 |
- * @endcode
- *
- * @ingroup LAOperators
- */
-
-// This is a workaround for a bug in <=gcc-4.7 that does not like partial
-// template default values in function definitions in combination with
-// local lambda expressions [1] in the function body. As a workaround
-// declare the function with all default types and parameters first such
-// that the function definition is without default types and parameters.
-//
-// [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
-
-template <unsigned int n,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
-lower_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &);
-
-
-template <unsigned int n, typename Range, typename Domain>
-LinearOperator<Range, Domain>
-lower_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &block_matrix)
-{
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_range_vector = [block_matrix](Range &v, bool fast)
-  {
-    v.reinit(n);
-
-    for (unsigned int i = 0; i < n; ++i)
-      block_matrix[i][i].reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.reinit_domain_vector = [block_matrix](Domain &v, bool fast)
-  {
-    v.reinit(n);
-
-    for (unsigned int i = 0; i < n; ++i)
-      block_matrix[i][i].reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.vmult = [block_matrix](Range &v, const Domain &u)
-  {
-    v.block(0) = 0;
-    for (unsigned int i = 1; i < n; ++i)
-      {
-        block_matrix[i][0].vmult(v.block(i), u.block(0));
-        for (unsigned int j = 1; j < i; ++j)
-          block_matrix[i][j].vmult_add(v.block(i), u.block(j));
-      }
-  };
-
-  return_op.vmult_add = [block_matrix](Range &v, const Domain &u)
-  {
-    for (unsigned int i = 1; i < n; ++i)
-      {
-        block_matrix[i][0].vmult_add(v.block(i), u.block(0));
-        for (unsigned int j = 1; j < i; ++j)
-          block_matrix[i][j].vmult_add(v.block(i), u.block(j));
-      }
-  };
-
-  return_op.Tvmult = [block_matrix](Domain &v, const Range &u)
-  {
-    for (unsigned int i = 0; i < n; ++i)
-      {
-        v.block(i) = 0;
-        for (unsigned int j = i+1; j < n; ++j)
-          block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
-      }
-  };
-
-  return_op.Tvmult_add = [block_matrix](Domain &v, const Range &u)
-  {
-    for (unsigned int i = 0; i < n; ++i)
-      for (unsigned int j = i+1; j < n; ++j)
-        block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
-  };
-
-  return return_op;
-}
-
-/**
- * @relates LinearOperator
- *
- * This function is a specification of the above function that
- * allows to work with block matrices @p block_matrix .
- *
- * @ingroup LAOperators
- */
-
-template <unsigned int n,
-          typename Range = BlockVector<double>,
-          typename Domain = Range,
-          typename BlockMatrix>
-LinearOperator<Range, Domain>
-lower_triangular_operator(const BlockMatrix &);
-
-
-template <unsigned int n, typename Range, typename Domain, typename BlockMatrix>
-LinearOperator<Range, Domain>
-lower_triangular_operator(const BlockMatrix &block_matrix)
-{
-  Assert(block_matrix.n_block_rows() == block_matrix.n_block_cols(),
-         ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) );
-
-  std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> M;
-  for (unsigned int i = 0; i<n; ++i)
-    {
-      for (unsigned int j = 0; j<i; ++j)
-        M[i][j] = linear_operator<typename Range::BlockType, typename Domain::BlockType>(block_matrix.block(i,j));
-      M[i][i] = null_operator(linear_operator<typename Range::BlockType, typename Domain::BlockType>(block_matrix.block(i,i)).reinit_range_vector);
-    }
-  return lower_triangular_operator<n, Range, Domain>(M);
-}
-
-/**
- * @relates LinearOperator
- *
- * A function that takes an array of arrays of LinearOperators @p block_matrix
- * and returns its associated upper triangular matrix operator
- * (diagonal is not included).
- *
- * @code
- *  a00 | a01 | a02                 | a01 | a02
- *  ---------------             ---------------
- *  a10 | a11 | a12     ->          |     | a12
- *  ---------------             ---------------
- *  a20 | a21 | a22                 |     |
- * @endcode
- *
- * @ingroup LAOperators
- */
-
-// This is a workaround for a bug in <=gcc-4.7 that does not like partial
-// template default values in function definitions in combination with
-// local lambda expressions [1] in the function body. As a workaround
-// declare the function with all default types and parameters first such
-// that the function definition is without default types and parameters.
-//
-// [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
-
-template <unsigned int n,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
-upper_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &);
-
-
-template <unsigned int n, typename Range, typename Domain>
-LinearOperator<Range, Domain>
-upper_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &block_matrix)
-{
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_range_vector = [block_matrix](Range &v, bool fast)
-  {
-    v.reinit(n);
-
-    for (unsigned int i = 0; i < n; ++i)
-      block_matrix[i][i].reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.reinit_domain_vector = [block_matrix](Domain &v, bool fast)
-  {
-    v.reinit(n);
-
-    for (unsigned int i = 0; i < n; ++i)
-      block_matrix[i][i].reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.vmult = [block_matrix](Range &v, const Domain &u)
-  {
-    for (unsigned int i = 0; i < n - 1; ++i)
-      {
-        block_matrix[i][n - 1].vmult(v.block(i), u.block(n - 1));
-        for (unsigned int j = n - 2; j > i; --j)
-          block_matrix[i][j].vmult_add(v.block(i), u.block(j));
-      }
-    v.block(n - 1) = 0;
-  };
-
-  return_op.vmult_add = [block_matrix](Range &v, const Domain &u)
-  {
-    for (unsigned int i = 0; i < n; ++i)
-      for (unsigned int j = i + 1; j < n; ++j)
-        block_matrix[i][j].vmult_add(v.block(i), u.block(j));
-  };
-
-  return_op.Tvmult = [block_matrix](Domain &v, const Range &u)
-  {
-    for (unsigned int i = 0; i < n; ++i)
-      {
-        v.block(i) = 0;
-        for (unsigned int j = 0; j < i; ++j)
-          block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
-      }
-  };
-
-  return_op.Tvmult_add = [block_matrix](Domain &v, const Range &u)
-  {
-    for (unsigned int i = 0; i < n; ++i)
-      for (unsigned int j = 0; j < i; ++j)
-        block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
-  };
-
-  return return_op;
-}
-
-
-/**
- * @relates LinearOperator
- *
- * This function is a specification of the above function that
- * allows to work with block matrices @p block_matrix .
- *
- * @ingroup LAOperators
- */
-
-template <unsigned int n,
-          typename Range = BlockVector<double>,
-          typename Domain = Range,
-          typename BlockMatrix>
-LinearOperator<Range, Domain>
-upper_triangular_operator(const BlockMatrix &);
-
-
-template <unsigned int n, typename Range, typename Domain, typename BlockMatrix>
-LinearOperator<Range, Domain>
-upper_triangular_operator(const BlockMatrix &block_matrix)
-{
-  Assert(block_matrix.n_block_rows() == block_matrix.n_block_cols(),
-         ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) );
-
-  std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> M;
-  for (unsigned int i = 0; i<n; ++i)
-    {
-      M[i][i] = null_operator(linear_operator<typename Range::BlockType, typename Domain::BlockType>(block_matrix.block(i,i)).reinit_range_vector);
-
-      for (unsigned int j = i+1; j<n; ++j)
-        M[i][j] = linear_operator<typename Range::BlockType, typename Domain::BlockType>(block_matrix.block(i,j));
-    }
-  return upper_triangular_operator<n, Range, Domain>(M);
-}
-
-
-/**
- * @relates LinearOperator
- *
- * Let M be a block matrix of the form Id + T made of nxn blocks and where
- * T is a lower / upper triangular (without diagonal). Then, its inverse is
- * of the form:
- * @code
- *  Id + sum_{i=1}^{n-1} (-1)^i T^i
- * @endcode
- * This formula can be used to invert all triangular matrices (diagonal
- * included of course).
- *
- * This function takes a block matrix @p block_matrix (possibly full) and a
- * linear block operator  @p inverse_diagonal made of the inverse of the
- * diagonal blocks inverses. The output is the inverse of the matrix in the
- * case of a triangular matrix and as inverse_diagonal its diagonal blocks
- * inverses. Otherwise, the result is a preconditioner.
- *
- * The parameter @p lower is a bool that allows to specify if we want to
- * use lower triangular part of @p block_matrix (true, this is the default
- * value) or to use upper triangular part of @p block_matrix (false).
- *
- * @ingroup LAOperators
- */
-
-// workaround for a bug in <=gcc-4.7 that does not like partial template
-// default values in combination with local lambda expressions [1]
-// [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
-template <unsigned int n,
-          typename Range = BlockVector<double>,
-          typename Domain = Range,
-          typename BlockMatrix>
-LinearOperator<Range, Domain>
-block_triangular_inverse(const BlockMatrix &,
-                         const LinearOperator<Range, Domain> &,
-                         bool lower = true);
-
-
-template <unsigned int n, typename Range, typename Domain, typename BlockMatrix>
-LinearOperator<Range, Domain>
-block_triangular_inverse(const BlockMatrix &block_matrix,
-                         const LinearOperator<Range, Domain> &inverse_diagonal,
-                         bool lower)
-{
-  Assert(block_matrix.n_block_rows() == n,
-         ExcDimensionMismatch(block_matrix.n_block_rows(), n));
-  Assert(block_matrix.n_block_rows() == block_matrix.n_block_cols(),
-         ExcDimensionMismatch(block_matrix.n_block_rows(),
-                              block_matrix.n_block_cols()));
-
-  LinearOperator<Range, Domain> op_a;
-
-  if (lower)
-    {
-      op_a = lower_triangular_operator<n, Range, Domain, BlockMatrix>(block_matrix);
-    }
-  else
-    {
-      op_a = upper_triangular_operator<n, Range, Domain, BlockMatrix>(block_matrix);
-    }
-
-  auto id     = identity_operator(op_a.reinit_range_vector);
-  auto result = identity_operator(op_a.reinit_range_vector);
-
-  // Notice that the following formula is recursive. We are evaluating:
-  // Id - T + T^2 - T^3 ... (- T)^n
-  for (unsigned int i = 0; i < n - 1; ++i)
-    result = id - inverse_diagonal * op_a * result;
-
-  return result * inverse_diagonal;
-}
-
 /**
  * @relates LinearOperator
  *
  * This function implement a forward substitution argument to invert a lower
  * block triangular matrix.
- * It takes as argement a block matrix @p block_matrix and an array of LinearOperators
+ * It takes as argement an array of array of LinearOperators @p block_matrix
+ * representing a block matrix and an array of LinearOperators
  * @p inverse_diagonal representing inverses of digonal blocks of @p block_matrix.
  *
  * Let us assume we have a linear system where each coefficient of the system is a
@@ -686,12 +347,16 @@ block_triangular_inverse(const BlockMatrix &block_matrix,
  * Caveat: Tvmult and Tvmult_add have not been implemented, yet. This may lead to mistakes.
  * @ingroup LAOperators
  */
+template < unsigned int n,
+           typename Range = BlockVector<double>>
+LinearOperator<Range, Range>
+block_forward_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &,
+                           const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &);
 
-template <typename BlockMatrix,
-          unsigned int n  = BlockMatrix::n ,
-          typename Range  = typename BlockMatrix::BlockType>
+template < unsigned int n,
+           typename Range>
 LinearOperator<Range, Range>
-block_forward_substitution(const BlockMatrix &block_matrix,
+block_forward_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &block_matrix,
                            const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal)
 {
   LinearOperator<Range, Range> return_op;
@@ -720,7 +385,7 @@ block_forward_substitution(const BlockMatrix &block_matrix,
     v.collect_sizes();
   };
 
-  return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Range &u)
+  return_op.vmult = [block_matrix, inverse_diagonal](Range &v, const Range &u)
   {
     static GrowingVectorMemory<typename  Range::BlockType> vector_memory;
     typename  Range::BlockType *tmp = vector_memory.alloc();
@@ -732,7 +397,7 @@ block_forward_substitution(const BlockMatrix &block_matrix,
         *tmp = u.block(i);
         *tmp *= -1.;
         for (unsigned int j=0; j<i; ++j)
-          block_matrix.block(i,j).vmult_add(*tmp, v.block(j));
+          block_matrix[i][j].vmult_add(*tmp, v.block(j));
         *tmp *= -1.;
         inverse_diagonal[i].vmult(v.block(i),*tmp);
       }
@@ -740,7 +405,7 @@ block_forward_substitution(const BlockMatrix &block_matrix,
     vector_memory.free(tmp);
   };
 
-  return_op.vmult_add = [&block_matrix, inverse_diagonal](Range &v, const Range &u)
+  return_op.vmult_add = [block_matrix, inverse_diagonal](Range &v, const Range &u)
   {
     static GrowingVectorMemory<typename  Range::BlockType> vector_memory;
     typename  Range::BlockType *tmp = vector_memory.alloc();
@@ -752,7 +417,7 @@ block_forward_substitution(const BlockMatrix &block_matrix,
         *tmp = u.block(i);
         *tmp *= -1.;
         for (unsigned int j=0; j<i; ++j)
-          block_matrix.block(i,j).vmult_add(*tmp, v.block(j));
+          block_matrix[i][j].vmult_add(*tmp, v.block(j));
         *tmp *= -1.;
         inverse_diagonal[i].vmult_add(v.block(i),*tmp);
       }
@@ -768,7 +433,8 @@ block_forward_substitution(const BlockMatrix &block_matrix,
  *
  * This function implement a back substitution argument to invert an upper
  * block triangular matrix.
- * It takes as argement a block matrix @p block_matrix and an array of LinearOperators
+ * It takes as argement an array of array of LinearOperators @p block_matrix
+ * representing a block matrix and an array of LinearOperators
  * @p inverse_diagonal representing inverses of digonal blocks of @p block_matrix.
  *
  * Let us assume we have a linear system where each coefficient of the system is a
@@ -788,12 +454,16 @@ block_forward_substitution(const BlockMatrix &block_matrix,
  * Caveat: Tvmult and Tvmult_add have not been implemented, yet. This may lead to mistakes.
  * @ingroup LAOperators
  */
+template < unsigned int n,
+           typename Range = BlockVector<double>>
+LinearOperator<Range, Range>
+block_back_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &,
+                        const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &);
 
-template <typename BlockMatrix,
-          unsigned int n  = BlockMatrix::n ,
-          typename Range  = typename BlockMatrix::BlockType>
+template < unsigned int n,
+           typename Range>
 LinearOperator<Range, Range>
-block_back_substitution(const BlockMatrix &block_matrix,
+block_back_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &block_matrix,
                         const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal)
 {
   LinearOperator<Range, Range> return_op;
@@ -822,7 +492,7 @@ block_back_substitution(const BlockMatrix &block_matrix,
     v.collect_sizes();
   };
 
-  return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Range &u)
+  return_op.vmult = [block_matrix, inverse_diagonal](Range &v, const Range &u)
   {
     static GrowingVectorMemory<typename  Range::BlockType> vector_memory;
     typename  Range::BlockType *tmp = vector_memory.alloc();
@@ -834,7 +504,7 @@ block_back_substitution(const BlockMatrix &block_matrix,
         *tmp = u.block(i);
         *tmp *= -1.;
         for (int j=i+1; j<n; ++j)
-          block_matrix.block(i,j).vmult_add(*tmp,v.block(j));
+          block_matrix[i][j].vmult_add(*tmp,v.block(j));
         *tmp *= -1.;
         inverse_diagonal[i].vmult(v.block(i),*tmp);
       }
@@ -854,7 +524,7 @@ block_back_substitution(const BlockMatrix &block_matrix,
         *tmp = u.block(i);
         *tmp *= -1.;
         for (int j=i+1; j<n; ++j)
-          block_matrix.block(i,j).vmult_add(*tmp,v.block(j));
+          block_matrix[i][j].vmult_add(*tmp,v.block(j));
         *tmp *= -1.;
         inverse_diagonal[i].vmult_add(v.block(i),*tmp);
       }
@@ -865,6 +535,55 @@ block_back_substitution(const BlockMatrix &block_matrix,
   return return_op;
 }
 
+/**
+ * @relates LinearOperator
+ *
+ * This function uses above functions block_back_substitution and block_forward_substitution
+ * to invert triangular matrices.
+ * It takes as input a triangular block matrix @p block_matrix, an array of LinearOperators
+ * @p inverse_diagonal representing inverses of block_matrix, and an optional bool @p lower
+ * used to specify if block_matrix should be conidered as lower triangular matrix (true) or
+ * as upper triangular matrix (false). @p lower is equal to true by default.
+ *
+ */
+
+// workaround for a bug in <=gcc-4.7 that does not like partial template
+// default values in combination with local lambda expressions [1]
+// [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
+
+template <unsigned int n,
+          typename BlockMatrix,
+          typename Range = BlockVector<double>>
+LinearOperator<Range, Range>
+block_triangular_inverse(const BlockMatrix &,
+                         const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &,
+                         bool lower = true);
+
+template <unsigned int n,
+          typename BlockMatrix,
+          typename Range>
+LinearOperator<Range, Range>
+block_triangular_inverse(const BlockMatrix &block_matrix,
+                         const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal,
+                         bool lower)
+{
+  Assert(block_matrix.n_block_rows() == n,
+         ExcDimensionMismatch(block_matrix.n_block_rows(), n));
+  Assert(block_matrix.n_block_rows() == block_matrix.n_block_cols(),
+         ExcDimensionMismatch(block_matrix.n_block_rows(),
+                              block_matrix.n_block_cols()));
+
+  std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> M;
+  for (unsigned int i = 0; i<n; ++i)
+    for (unsigned int j = 0; j<n; ++j)
+      M[i][j] = linear_operator<typename Range::BlockType, typename Range::BlockType>(block_matrix.block(i,j));
+
+  if (lower)
+    return block_forward_substitution<n, Range>(M, inverse_diagonal);
+  else
+    return block_back_substitution<n, Range>(M, inverse_diagonal);
+}
+
 //@}
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/lac/block_linear_operator_03.cc b/tests/lac/block_linear_operator_03.cc
new file mode 100644 (file)
index 0000000..dad87eb
--- /dev/null
@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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 lower_triangular_operator:
+
+#include "../tests.h"
+
+#include <deal.II/lac/block_linear_operator.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#define PRINTME(title, name, var) \
+  deallog << title << std::endl; \
+  deallog << "Block vector: " name << ":" << std::endl; \
+  for (unsigned int i = 0; i < var.n_blocks(); ++i) \
+    deallog << "[block " << i << " ]  " << var.block(i);
+
+
+using namespace dealii;
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(2);
+
+  // BlockSparseMatrix:
+  {
+    BlockDynamicSparsityPattern dsp(3, 3);
+    for (unsigned int i = 0; i < 3; ++i)
+      for (unsigned int j = 0; j < 3; ++j)
+        dsp.block(i, j).reinit (1, 1);
+    dsp.collect_sizes ();
+
+    BlockSparsityPattern sparsity_pattern;
+    sparsity_pattern.copy_from(dsp);
+    sparsity_pattern.compress();
+
+    //  | 2 1 0 |     | 0 0 0 |
+    //  | 3 2 1 | ->  | 3 0 0 |
+    //  | 4 3 2 |     | 4 3 0 |
+    BlockSparseMatrix<double> a (sparsity_pattern);
+    for (unsigned int i = 0; i < 3; ++i)
+      for (unsigned int j = 0; j < 3; ++j)
+        a.block(i,j).set(0, 0, 2 + i - j);
+
+    auto op_a = linear_operator<BlockVector<double>>(a);
+    auto triangular_block_op = lower_triangular_operator<3, BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a);
+
+    BlockVector<double> u;
+    op_a.reinit_domain_vector(u, false);
+    BlockVector<double> v;
+    op_a.reinit_range_vector(v, false);
+
+
+    // vmult:
+    for(unsigned int i = 0; i<3; ++i)
+    {
+      u.block(i)[0] = i+1;
+      v.block(i)[0] = 1;
+    }
+
+    triangular_block_op.vmult(v, u);
+    PRINTME(" -- vmult --", "v", v);
+
+    // vmult_add:
+    for(unsigned int i = 0; i<3; ++i)
+      v.block(i)[0] = 1;
+
+    triangular_block_op.vmult_add(v, u);
+    PRINTME(" -- vmult_add --", "v", v);
+
+    // Tvmult
+    for(unsigned int i = 0; i<3; ++i)
+      v.block(i)[0] = i+1;
+
+    triangular_block_op.Tvmult(u, v);
+    PRINTME(" -- Tvmult --", "u", u);
+
+    // Tvmult_add
+    for(unsigned int i = 0; i<3; ++i)
+      u.block(i)[0] = 1;
+
+    triangular_block_op.Tvmult_add(u, v);
+    PRINTME(" -- Tvmult_add --", "u", u);
+  }
+}
similarity index 96%
rename from tests/lac/block_linear_operator_09.cc
rename to tests/lac/block_linear_operator_04.cc
index ad4532af790f586ff1ae7df5ec56ea81f2c571c5..e49548d5089ea0ce0a4e285fcf97ffdd466fa2d4 100644 (file)
@@ -66,7 +66,7 @@ int main()
     auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
     auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
 
-    auto inverse_op_a = block_forward_substitution< BlockSparseMatrix<double> >(a, {{d00, d11, d22}});
+    auto inverse_op_a = block_triangular_inverse<3, BlockSparseMatrix<double > >(a, {{d00, d11, d22}});
     auto identity = inverse_op_a * op_a;
 
     BlockVector<double> u;
@@ -153,7 +153,7 @@ int main()
     auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
     auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
 
-    auto inverse_op_a = block_back_substitution< BlockSparseMatrix<double> >(a, {{d00, d11, d22}});
+    auto inverse_op_a = block_triangular_inverse< 3, BlockSparseMatrix<double > >(a, {{d00, d11, d22}}, false);
     auto identity = inverse_op_a * op_a;
 
     BlockVector<double> u;
diff --git a/tests/lac/block_linear_operator_05.cc b/tests/lac/block_linear_operator_05.cc
new file mode 100644 (file)
index 0000000..7c8378c
--- /dev/null
@@ -0,0 +1,249 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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 block_back_substitution and block_forward_substitution:
+
+#include "../tests.h"
+
+#include <deal.II/lac/block_linear_operator.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#define PRINTME(name, var) \
+  deallog << "Block vector: " name << ":" << std::endl; \
+  for (unsigned int i = 0; i < var.n_blocks(); ++i) \
+    deallog << "[block " << i << " ]  " << var.block(i);
+
+
+using namespace dealii;
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(12);
+
+  // BlockSparseMatrix:
+  {
+    BlockDynamicSparsityPattern dsp(3, 3);
+    for (unsigned int i = 0; i < 3; ++i)
+      for (unsigned int j = 0; j < 3; ++j)
+        dsp.block(i, j).reinit (1, 1);
+    dsp.collect_sizes ();
+
+    BlockSparsityPattern sparsity_pattern;
+    sparsity_pattern.copy_from(dsp);
+    sparsity_pattern.compress();
+
+    BlockSparseMatrix<double> a (sparsity_pattern);
+    for (unsigned int i = 0; i < 3; ++i)
+    {
+      a.block(i,i).set(0, 0, i+i +1);
+      for (unsigned int j = 0; j < i; ++j)
+        a.block(i,j).set(0, 0, 10);
+    }
+
+    BlockSparseMatrix<double> d(sparsity_pattern);
+    for (unsigned int i = 0; i < 3; ++i)
+        d.block(i,i).set(0,0, 1.0 / (i+i +1) );
+
+    auto op_a         = linear_operator< BlockVector<double> >(a);
+
+    auto a00 = linear_operator< Vector<double>, Vector<double> >(a.block(0,0));
+    auto a01 = linear_operator< Vector<double>, Vector<double> >(a.block(0,1));
+    auto a02 = linear_operator< Vector<double>, Vector<double> >(a.block(0,2));
+    auto a10 = linear_operator< Vector<double>, Vector<double> >(a.block(1,0));
+    auto a11 = linear_operator< Vector<double>, Vector<double> >(a.block(1,1));
+    auto a12 = linear_operator< Vector<double>, Vector<double> >(a.block(1,2));
+    auto a20 = linear_operator< Vector<double>, Vector<double> >(a.block(2,0));
+    auto a21 = linear_operator< Vector<double>, Vector<double> >(a.block(2,1));
+    auto a22 = linear_operator< Vector<double>, Vector<double> >(a.block(2,2));
+
+    auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
+    auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
+    auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
+
+    auto inverse_op_a = block_forward_substitution< 3, BlockVector<double> >(
+      {{
+        {{a00, a01, a02}},
+        {{a10, a11, a12}},
+        {{a20, a21, a22}}
+      }},
+      { {d00, d11, d22}});
+
+    auto identity = inverse_op_a * op_a;
+
+    BlockVector<double> u;
+    BlockVector<double> v;
+
+    deallog << " -- Matrix -- " << std::endl;
+    op_a.reinit_domain_vector(u, false);
+    op_a.reinit_range_vector(v, false);
+    for(unsigned int j = 0; j<3; ++j)
+    {
+      for(unsigned int i = 0; i<3; ++i)
+      {
+        u.block(i)[0] = 0;;
+        v.block(i)[0] = 0;
+      }
+      u.block(j)[0] = 1;
+
+      op_a.vmult(v, u);
+
+      PRINTME("v", v);
+    }
+
+    deallog << " -- Inverse -- " << std::endl;
+    inverse_op_a.reinit_domain_vector(u, false);
+    inverse_op_a.reinit_range_vector(v, true);
+    for(unsigned int j = 0; j<3; ++j)
+    {
+      for(unsigned int i = 0; i<3; ++i)
+      {
+        u.block(i)[0] = 0;
+        v.block(i)[0] = 0;
+      }
+      u.block(j)[0] = 1;;
+
+      inverse_op_a.vmult(v, u);
+
+      PRINTME("v", v);
+    }
+
+    deallog << " -- Identity -- " << std::endl;
+    identity.reinit_domain_vector(u, false);
+    identity.reinit_range_vector(v, false);
+    for(unsigned int j = 0; j<3; ++j)
+    {
+      for(unsigned int i = 0; i<3; ++i)
+      {
+        u.block(i)[0] = 0;;
+        v.block(i)[0] = 0;
+      }
+      u.block(j)[0] = 1;;
+
+      identity.vmult(v, u);
+
+      PRINTME("v", v);
+    }
+  }
+
+
+  {
+    BlockDynamicSparsityPattern dsp(3, 3);
+    for (unsigned int i = 0; i < 3; ++i)
+      for (unsigned int j = 0; j < 3; ++j)
+        dsp.block(i, j).reinit (1, 1);
+    dsp.collect_sizes ();
+
+    BlockSparsityPattern sparsity_pattern;
+    sparsity_pattern.copy_from(dsp);
+    sparsity_pattern.compress();
+
+    BlockSparseMatrix<double> a (sparsity_pattern);
+    for (unsigned int i = 0; i < 3; ++i)
+    {
+      a.block(i,i).set(0, 0, i+i +1);
+      for (unsigned int j = i+1; j < 3; ++j)
+        a.block(i,j).set(0, 0, 10);
+    }
+
+    BlockSparseMatrix<double> d(sparsity_pattern);
+    for (unsigned int i = 0; i < 3; ++i)
+        d.block(i,i).set(0,0, 1.0 / (i+i +1) );
+
+        auto op_a         = linear_operator< BlockVector<double> >(a);
+
+        auto a00 = linear_operator< Vector<double>, Vector<double> >(a.block(0,0));
+        auto a01 = linear_operator< Vector<double>, Vector<double> >(a.block(0,1));
+        auto a02 = linear_operator< Vector<double>, Vector<double> >(a.block(0,2));
+        auto a10 = linear_operator< Vector<double>, Vector<double> >(a.block(1,0));
+        auto a11 = linear_operator< Vector<double>, Vector<double> >(a.block(1,1));
+        auto a12 = linear_operator< Vector<double>, Vector<double> >(a.block(1,2));
+        auto a20 = linear_operator< Vector<double>, Vector<double> >(a.block(2,0));
+        auto a21 = linear_operator< Vector<double>, Vector<double> >(a.block(2,1));
+        auto a22 = linear_operator< Vector<double>, Vector<double> >(a.block(2,2));
+
+        auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
+        auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
+        auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
+
+        auto inverse_op_a = block_back_substitution< 3, BlockVector<double> >(
+          {{
+            {{a00, a01, a02}},
+            {{a10, a11, a12}},
+            {{a20, a21, a22}}
+          }},
+          { {d00, d11, d22}});
+
+        auto identity = inverse_op_a * op_a;
+
+    BlockVector<double> u;
+    BlockVector<double> v;
+
+    deallog << " -- Matrix -- " << std::endl;
+    op_a.reinit_domain_vector(u, false);
+    op_a.reinit_range_vector(v, false);
+    for(unsigned int j = 0; j<3; ++j)
+    {
+      for(unsigned int i = 0; i<3; ++i)
+      {
+        u.block(i)[0] = 0;;
+        v.block(i)[0] = 0;
+      }
+      u.block(j)[0] = 1;
+
+      op_a.vmult(v, u);
+
+      PRINTME("v", v);
+    }
+
+    deallog << " -- Inverse -- " << std::endl;
+    inverse_op_a.reinit_domain_vector(u, false);
+    inverse_op_a.reinit_range_vector(v, true);
+    for(unsigned int j = 0; j<3; ++j)
+    {
+      for(unsigned int i = 0; i<3; ++i)
+      {
+        u.block(i)[0] = 0;
+        v.block(i)[0] = 0;
+      }
+      u.block(j)[0] = 1;;
+
+      inverse_op_a.vmult(v, u);
+
+      PRINTME("v", v);
+    }
+
+    deallog << " -- Identity -- " << std::endl;
+    identity.reinit_domain_vector(u, false);
+    identity.reinit_range_vector(v, false);
+    for(unsigned int j = 0; j<3; ++j)
+    {
+      for(unsigned int i = 0; i<3; ++i)
+      {
+        u.block(i)[0] = 0;;
+        v.block(i)[0] = 0;
+      }
+      u.block(j)[0] = 1;;
+
+      identity.vmult(v, u);
+
+      PRINTME("v", v);
+    }
+  }
+}
diff --git a/tests/lac/block_linear_operator_05.with_cxx11=on.output b/tests/lac/block_linear_operator_05.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..1e80e98
--- /dev/null
@@ -0,0 +1,79 @@
+
+DEAL:: -- Matrix -- 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  1.00000000000 
+DEAL::[block 1 ]  10.0000000000 
+DEAL::[block 2 ]  10.0000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  3.00000000000 
+DEAL::[block 2 ]  10.0000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  5.00000000000 
+DEAL:: -- Inverse -- 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  1.00000000000 
+DEAL::[block 1 ]  -3.33333333333 
+DEAL::[block 2 ]  4.66666666667 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  0.333333333333 
+DEAL::[block 2 ]  -0.666666666667 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  0.200000000000 
+DEAL:: -- Identity -- 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  1.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  1.00000000000 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  1.00000000000 
+DEAL:: -- Matrix -- 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  1.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  10.0000000000 
+DEAL::[block 1 ]  3.00000000000 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  10.0000000000 
+DEAL::[block 1 ]  10.0000000000 
+DEAL::[block 2 ]  5.00000000000 
+DEAL:: -- Inverse -- 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  1.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  -3.33333333333 
+DEAL::[block 1 ]  0.333333333333 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  4.66666666667 
+DEAL::[block 1 ]  -0.666666666667 
+DEAL::[block 2 ]  0.200000000000 
+DEAL:: -- Identity -- 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  1.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  1.00000000000 
+DEAL::[block 2 ]  0.00000000000 
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.00000000000 
+DEAL::[block 1 ]  0.00000000000 
+DEAL::[block 2 ]  1.00000000000 
diff --git a/tests/lac/block_linear_operator_06.cc b/tests/lac/block_linear_operator_06.cc
deleted file mode 100644 (file)
index 69e9b65..0000000
+++ /dev/null
@@ -1,125 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2015 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 preconditioner_from_diagonal_inverse:
-
-#include "../tests.h"
-
-#include <deal.II/lac/block_linear_operator.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#define PRINTME(name, var) \
-  deallog << "Block vector: " name << ":" << std::endl; \
-  for (unsigned int i = 0; i < var.n_blocks(); ++i) \
-    deallog << "[block " << i << " ]  " << var.block(i);
-
-
-using namespace dealii;
-
-int main()
-{
-  initlog();
-  deallog << std::setprecision(12);
-
-  // BlockSparseMatrix:
-  {
-    BlockDynamicSparsityPattern dsp(3, 3);
-    for (unsigned int i = 0; i < 3; ++i)
-      for (unsigned int j = 0; j < 3; ++j)
-        dsp.block(i, j).reinit (1, 1);
-    dsp.collect_sizes ();
-
-    BlockSparsityPattern sparsity_pattern;
-    sparsity_pattern.copy_from(dsp);
-    sparsity_pattern.compress();
-
-    BlockSparseMatrix<double> a (sparsity_pattern);
-    for (unsigned int i = 0; i < 3; ++i)
-    {
-      a.block(i,i).set(0, 0, i+i +1);
-      for (unsigned int j = 0; j < i; ++j)
-        a.block(i,j).set(0, 0, 10);
-    }
-
-    BlockSparseMatrix<double> d(sparsity_pattern);
-    for (unsigned int i = 0; i < 3; ++i)
-        d.set(i,i, 1.0 / (i+i +1) );
-
-    auto op_a         = linear_operator< BlockVector<double> >(a);
-    auto diagonal_inv = linear_operator< BlockVector<double>, BlockVector<double> >(d);
-    // auto diagonal_inv = block_diagonal_operator(diag);
-    auto inverse_op_a = block_triangular_inverse<3, BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a, diagonal_inv);
-
-    auto identity = inverse_op_a * op_a;
-
-    BlockVector<double> u;
-    BlockVector<double> v;
-
-    deallog << " -- Matrix -- " << std::endl;
-    op_a.reinit_domain_vector(u, false);
-    op_a.reinit_range_vector(v, false);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;
-
-      op_a.vmult(v, u);
-
-      PRINTME("v", v);
-    }
-
-    deallog << " -- Inverse -- " << std::endl;
-    inverse_op_a.reinit_domain_vector(u, false);
-    inverse_op_a.reinit_range_vector(v, true);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;;
-
-      inverse_op_a.vmult(v, u);
-
-      PRINTME("v", v);
-    }
-
-    deallog << " -- Identity -- " << std::endl;
-    identity.reinit_domain_vector(u, false);
-    identity.reinit_range_vector(v, false);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;;
-
-      identity.vmult(v, u);
-
-      PRINTME("v", v);
-    }
-  }
-}

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.