From 0d50e5fe9259f16ed72b14510857e5c959cf1951 Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Tue, 16 Jun 2015 21:57:23 +0200 Subject: [PATCH] removing file depending on upper and lower triangular --- include/deal.II/lac/block_linear_operator.h | 433 +++--------------- tests/lac/block_linear_operator_03.cc | 101 ++++ ...k_linear_operator_03.with_cxx11=on.output} | 0 ...ator_09.cc => block_linear_operator_04.cc} | 4 +- ...k_linear_operator_04.with_cxx11=on.output} | 0 tests/lac/block_linear_operator_05.cc | 249 ++++++++++ ...ck_linear_operator_05.with_cxx11=on.output | 79 ++++ tests/lac/block_linear_operator_06.cc | 125 ----- 8 files changed, 507 insertions(+), 484 deletions(-) create mode 100644 tests/lac/block_linear_operator_03.cc rename tests/lac/{block_linear_operator_06.with_cxx11=on.output => block_linear_operator_03.with_cxx11=on.output} (100%) rename tests/lac/{block_linear_operator_09.cc => block_linear_operator_04.cc} (96%) rename tests/lac/{block_linear_operator_09.with_cxx11=on.output => block_linear_operator_04.with_cxx11=on.output} (100%) create mode 100644 tests/lac/block_linear_operator_05.cc create mode 100644 tests/lac/block_linear_operator_05.with_cxx11=on.output delete mode 100644 tests/lac/block_linear_operator_06.cc diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h index 1a480a7d45..64ff277b8f 100644 --- a/include/deal.II/lac/block_linear_operator.h +++ b/include/deal.II/lac/block_linear_operator.h @@ -322,352 +322,13 @@ block_diagonal_operator(const LinearOperator 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 , - typename Domain = Range> -LinearOperator -lower_triangular_operator(const std::array, n>, n> &); - - -template -LinearOperator -lower_triangular_operator(const std::array, n>, n> &block_matrix) -{ - LinearOperator 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 , - typename Domain = Range, - typename BlockMatrix> -LinearOperator -lower_triangular_operator(const BlockMatrix &); - - -template -LinearOperator -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, n>, n> M; - for (unsigned int i = 0; i(block_matrix.block(i,j)); - M[i][i] = null_operator(linear_operator(block_matrix.block(i,i)).reinit_range_vector); - } - return lower_triangular_operator(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 , - typename Domain = Range> -LinearOperator -upper_triangular_operator(const std::array, n>, n> &); - - -template -LinearOperator -upper_triangular_operator(const std::array, n>, n> &block_matrix) -{ - LinearOperator 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 , - typename Domain = Range, - typename BlockMatrix> -LinearOperator -upper_triangular_operator(const BlockMatrix &); - - -template -LinearOperator -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, n>, n> M; - for (unsigned int i = 0; i(block_matrix.block(i,i)).reinit_range_vector); - - for (unsigned int j = i+1; j(block_matrix.block(i,j)); - } - return upper_triangular_operator(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 , - typename Domain = Range, - typename BlockMatrix> -LinearOperator -block_triangular_inverse(const BlockMatrix &, - const LinearOperator &, - bool lower = true); - - -template -LinearOperator -block_triangular_inverse(const BlockMatrix &block_matrix, - const LinearOperator &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 op_a; - - if (lower) - { - op_a = lower_triangular_operator(block_matrix); - } - else - { - op_a = upper_triangular_operator(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> +LinearOperator +block_forward_substitution(const std::array, n>, n> &, + const std::array, n> &); -template +template < unsigned int n, + typename Range> LinearOperator -block_forward_substitution(const BlockMatrix &block_matrix, +block_forward_substitution(const std::array, n>, n> &block_matrix, const std::array, n> &inverse_diagonal) { LinearOperator 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 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 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> +LinearOperator +block_back_substitution(const std::array, n>, n> &, + const std::array, n> &); -template +template < unsigned int n, + typename Range> LinearOperator -block_back_substitution(const BlockMatrix &block_matrix, +block_back_substitution(const std::array, n>, n> &block_matrix, const std::array, n> &inverse_diagonal) { LinearOperator 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 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> +LinearOperator +block_triangular_inverse(const BlockMatrix &, + const std::array, n> &, + bool lower = true); + +template +LinearOperator +block_triangular_inverse(const BlockMatrix &block_matrix, + const std::array, 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, n>, n> M; + for (unsigned int i = 0; i(block_matrix.block(i,j)); + + if (lower) + return block_forward_substitution(M, inverse_diagonal); + else + return block_back_substitution(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 index 0000000000..dad87eb770 --- /dev/null +++ b/tests/lac/block_linear_operator_03.cc @@ -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 +#include +#include +#include +#include +#include + +#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 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>(a); + auto triangular_block_op = lower_triangular_operator<3, BlockVector, BlockVector, BlockSparseMatrix >(a); + + BlockVector u; + op_a.reinit_domain_vector(u, false); + BlockVector 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); + } +} diff --git a/tests/lac/block_linear_operator_06.with_cxx11=on.output b/tests/lac/block_linear_operator_03.with_cxx11=on.output similarity index 100% rename from tests/lac/block_linear_operator_06.with_cxx11=on.output rename to tests/lac/block_linear_operator_03.with_cxx11=on.output diff --git a/tests/lac/block_linear_operator_09.cc b/tests/lac/block_linear_operator_04.cc similarity index 96% rename from tests/lac/block_linear_operator_09.cc rename to tests/lac/block_linear_operator_04.cc index ad4532af79..e49548d508 100644 --- a/tests/lac/block_linear_operator_09.cc +++ b/tests/lac/block_linear_operator_04.cc @@ -66,7 +66,7 @@ int main() auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); - auto inverse_op_a = block_forward_substitution< BlockSparseMatrix >(a, {{d00, d11, d22}}); + auto inverse_op_a = block_triangular_inverse<3, BlockSparseMatrix >(a, {{d00, d11, d22}}); auto identity = inverse_op_a * op_a; BlockVector u; @@ -153,7 +153,7 @@ int main() auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); - auto inverse_op_a = block_back_substitution< BlockSparseMatrix >(a, {{d00, d11, d22}}); + auto inverse_op_a = block_triangular_inverse< 3, BlockSparseMatrix >(a, {{d00, d11, d22}}, false); auto identity = inverse_op_a * op_a; BlockVector u; diff --git a/tests/lac/block_linear_operator_09.with_cxx11=on.output b/tests/lac/block_linear_operator_04.with_cxx11=on.output similarity index 100% rename from tests/lac/block_linear_operator_09.with_cxx11=on.output rename to tests/lac/block_linear_operator_04.with_cxx11=on.output diff --git a/tests/lac/block_linear_operator_05.cc b/tests/lac/block_linear_operator_05.cc new file mode 100644 index 0000000000..7c8378cb2a --- /dev/null +++ b/tests/lac/block_linear_operator_05.cc @@ -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 +#include +#include +#include +#include +#include + +#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 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 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 >(a); + + auto a00 = linear_operator< Vector, Vector >(a.block(0,0)); + auto a01 = linear_operator< Vector, Vector >(a.block(0,1)); + auto a02 = linear_operator< Vector, Vector >(a.block(0,2)); + auto a10 = linear_operator< Vector, Vector >(a.block(1,0)); + auto a11 = linear_operator< Vector, Vector >(a.block(1,1)); + auto a12 = linear_operator< Vector, Vector >(a.block(1,2)); + auto a20 = linear_operator< Vector, Vector >(a.block(2,0)); + auto a21 = linear_operator< Vector, Vector >(a.block(2,1)); + auto a22 = linear_operator< Vector, Vector >(a.block(2,2)); + + auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); + auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); + auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); + + auto inverse_op_a = block_forward_substitution< 3, BlockVector >( + {{ + {{a00, a01, a02}}, + {{a10, a11, a12}}, + {{a20, a21, a22}} + }}, + { {d00, d11, d22}}); + + auto identity = inverse_op_a * op_a; + + BlockVector u; + BlockVector 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 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 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 >(a); + + auto a00 = linear_operator< Vector, Vector >(a.block(0,0)); + auto a01 = linear_operator< Vector, Vector >(a.block(0,1)); + auto a02 = linear_operator< Vector, Vector >(a.block(0,2)); + auto a10 = linear_operator< Vector, Vector >(a.block(1,0)); + auto a11 = linear_operator< Vector, Vector >(a.block(1,1)); + auto a12 = linear_operator< Vector, Vector >(a.block(1,2)); + auto a20 = linear_operator< Vector, Vector >(a.block(2,0)); + auto a21 = linear_operator< Vector, Vector >(a.block(2,1)); + auto a22 = linear_operator< Vector, Vector >(a.block(2,2)); + + auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); + auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); + auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); + + auto inverse_op_a = block_back_substitution< 3, BlockVector >( + {{ + {{a00, a01, a02}}, + {{a10, a11, a12}}, + {{a20, a21, a22}} + }}, + { {d00, d11, d22}}); + + auto identity = inverse_op_a * op_a; + + BlockVector u; + BlockVector 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 index 0000000000..1e80e98cd0 --- /dev/null +++ b/tests/lac/block_linear_operator_05.with_cxx11=on.output @@ -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 index 69e9b65d5d..0000000000 --- a/tests/lac/block_linear_operator_06.cc +++ /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 -#include -#include -#include -#include -#include - -#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 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 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 >(a); - auto diagonal_inv = linear_operator< BlockVector, BlockVector >(d); - // auto diagonal_inv = block_diagonal_operator(diag); - auto inverse_op_a = block_triangular_inverse<3, BlockVector, BlockVector, BlockSparseMatrix >(a, diagonal_inv); - - auto identity = inverse_op_a * op_a; - - BlockVector u; - BlockVector 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); - } - } -} -- 2.39.5