From: ESeNonFossiIo Date: Wed, 20 May 2015 13:36:44 +0000 (+0200) Subject: upper triangular block matrix X-Git-Tag: v8.3.0-rc1~155^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=caaa7800d3b114a80748f8f3c4774482288872f7;p=dealii.git upper triangular block matrix upper triangular block matrix tests comments crossed tests small typos --- diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index 0f39f75138..c46a0a44e5 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -2229,6 +2229,106 @@ lower_triangular_operator(BlockMatrix &a) } +template , + typename Domain = Range, + typename BlockMatrix> +LinearOperator +upper_triangular_operator(BlockMatrix &a) +{ + Assert( a.n_block_rows() == a.n_block_cols(), + ExcDimensionMismatch(a.n_block_rows(),a.n_block_cols()) ); + + LinearOperator return_op; + + return_op.reinit_range_vector = [&a](Range &v, bool fast) + { + // Reinitialize the block vector to have the number of blocks + // equal to the number of row blocks of the matrix a. + v.reinit(a.n_block_rows()); + + // And reinitialize every individual block with reinit_range_vectors: + for (unsigned int i = 0; i < a.n_block_rows(); ++i) + linear_operator(a.block(i,0)).reinit_range_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.reinit_domain_vector = [&a](Domain &v, bool fast) + { + // Reinitialize the block vector to have the number of blocks + // equal to the number of coloumn blocks of the matrix a. + v.reinit(a.n_block_cols()); + + // And reinitialize every individual block with reinit_domain_vectors: + for (unsigned int i = 0; i < a.n_block_cols(); ++i) + linear_operator(a.block(0,i)).reinit_domain_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.vmult = [&a](Range &v, const Domain &u) + { + Assert( v.n_blocks() == a.n_block_rows(), + ExcDimensionMismatch(v.n_blocks(), a.n_block_rows())); + Assert( u.n_blocks() == a.n_block_cols(), + ExcDimensionMismatch(u.n_blocks(), a.n_block_cols())); + + for (unsigned int i = 0; i < a.n_block_rows(); ++i) + { + v.block(i) *= 0; + for (unsigned int j = i + 1; j < a.n_block_cols(); ++j) + a.block(i,j).Tvmult_add(v.block(i), u.block(j)); + } + + }; + + return_op.vmult_add = [&a](Range &v, const Domain &u) + { + Assert( v.n_blocks() == a.n_block_rows(), + ExcDimensionMismatch(v.n_blocks(), a.n_block_rows())); + Assert( u.n_blocks() == a.n_block_cols(), + ExcDimensionMismatch(u.n_blocks(), a.n_block_cols())); + + for (unsigned int i = 0; i < a.n_block_cols(); ++i) + for (unsigned int j = i + 1; j < a.n_block_rows(); ++j) + a.block(i,j).Tvmult_add(v.block(i), u.block(j)); + }; + + + + return_op.Tvmult = [&a](Domain &v, const Range &u) + { + Assert( v.n_blocks() == a.n_block_cols(), + ExcDimensionMismatch(v.n_blocks(), a.n_block_cols())); + Assert( u.n_blocks() == a.n_block_rows(), + ExcDimensionMismatch(u.n_blocks(), a.n_block_rows())); + + + + for (unsigned int i = 0; i < a.n_block_cols(); ++i) + { + v.block(i) *= 0; + for (unsigned int j = 0; j < i; ++j) + a.block(j,i).vmult_add(v.block(i), u.block(j)); + } + }; + + + return_op.Tvmult_add = [&a](Domain &v, const Range &u) + { + Assert( v.n_blocks() == a.n_block_cols(), + ExcDimensionMismatch(v.n_blocks(), a.n_block_cols())); + Assert( u.n_blocks() == a.n_block_rows(), + ExcDimensionMismatch(u.n_blocks(), a.n_block_rows())); + + for (unsigned int i = 0; i < a.n_block_cols(); ++i) + for (unsigned int j = 0; j < i; ++j) + a.block(j,i).vmult_add(v.block(i), u.block(j)); + }; + + return return_op; +} + //@} DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/linear_operator_12.cc b/tests/lac/linear_operator_12.cc new file mode 100644 index 0000000000..f76d503b95 --- /dev/null +++ b/tests/lac/linear_operator_12.cc @@ -0,0 +1,105 @@ +// --------------------------------------------------------------------- +// +// 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 non symmetric variants: + +#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(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 3 4 | + // | 1 2 3 | + // | 0 1 2 | + 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 + j - i); + + auto op_a = linear_operator>(a); + + auto triangular_block_op = upper_triangular_operator< 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); + deallog << " -- vmult --" << std::endl; + PRINTME("v", v); + + // vmult_add: + for(unsigned int i = 0; i<3; ++i) + v.block(i)[0] = 1; + + triangular_block_op.vmult_add(v, u); + deallog << " -- vmult_add --" << std::endl; + PRINTME("v", v); + + // Tvmult + for(unsigned int i = 0; i<3; ++i) + v.block(i)[0] = i+1; + + triangular_block_op.Tvmult(u, v); + deallog << " -- Tvmult --" << std::endl; + PRINTME("u", u); + + // Tvmult_add + for(unsigned int i = 0; i<3; ++i) + u.block(i)[0] = 1; + + triangular_block_op.Tvmult_add(u, v); + deallog << " -- Tvmult_add --" << std::endl; + PRINTME("u", u); + } +} diff --git a/tests/lac/linear_operator_12.with_cxx11=on.output b/tests/lac/linear_operator_12.with_cxx11=on.output new file mode 100644 index 0000000000..822b43771b --- /dev/null +++ b/tests/lac/linear_operator_12.with_cxx11=on.output @@ -0,0 +1,21 @@ + +DEAL:: -- vmult -- +DEAL::Block vector: v: +DEAL::[block 0 ] 18. +DEAL::[block 1 ] 9.0 +DEAL::[block 2 ] 0.0 +DEAL:: -- vmult_add -- +DEAL::Block vector: v: +DEAL::[block 0 ] 19. +DEAL::[block 1 ] 10. +DEAL::[block 2 ] 1.0 +DEAL:: -- Tvmult -- +DEAL::Block vector: u: +DEAL::[block 0 ] 0.0 +DEAL::[block 1 ] 3.0 +DEAL::[block 2 ] 10. +DEAL:: -- Tvmult_add -- +DEAL::Block vector: u: +DEAL::[block 0 ] 1.0 +DEAL::[block 1 ] 4.0 +DEAL::[block 2 ] 11. diff --git a/tests/lac/linear_operator_13.cc b/tests/lac/linear_operator_13.cc new file mode 100644 index 0000000000..199bd702d3 --- /dev/null +++ b/tests/lac/linear_operator_13.cc @@ -0,0 +1,135 @@ +// --------------------------------------------------------------------- +// +// 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 non symmetric variants: + +#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(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(); + + 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 + j + i); + + auto op_a = linear_operator>(a); + + auto lower_triangular_block_op = lower_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); + auto upper_triangular_block_op = upper_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); + + + BlockVector u1; + op_a.reinit_domain_vector(u1, false); + BlockVector v1; + op_a.reinit_range_vector(v1, false); + BlockVector u2; + op_a.reinit_domain_vector(u2, false); + BlockVector v2; + op_a.reinit_range_vector(v2, false); + + + // check 1 + for(unsigned int i = 0; i<3; ++i) + { + u1.block(i)[0] = i+1; + u2.block(i)[0] = i+1; + v1.block(i)[0] = 1; + v2.block(i)[0] = 1; + } + + lower_triangular_block_op.vmult(v1, u1); + upper_triangular_block_op.Tvmult(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + + + // check 2 + for(unsigned int i = 0; i<3; ++i) + { + u1.block(i)[0] = i+1; + u2.block(i)[0] = i+1; + v1.block(i)[0] = 1; + v2.block(i)[0] = 1; + } + + lower_triangular_block_op.Tvmult(v1, u1); + upper_triangular_block_op.vmult(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + + // check 3 + for(unsigned int i = 0; i<3; ++i) + { + u1.block(i)[0] = i+1; + u2.block(i)[0] = i+1; + v1.block(i)[0] = 1; + v2.block(i)[0] = 1; + } + + lower_triangular_block_op.vmult_add(v1, u1); + upper_triangular_block_op.Tvmult_add(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + + + // check 4 + for(unsigned int i = 0; i<3; ++i) + { + u1.block(i)[0] = i+1; + u2.block(i)[0] = i+1; + v1.block(i)[0] = 1; + v2.block(i)[0] = 1; + } + + lower_triangular_block_op.Tvmult_add(v1, u1); + upper_triangular_block_op.vmult_add(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + } +} diff --git a/tests/lac/linear_operator_13.with_cxx11=on.output b/tests/lac/linear_operator_13.with_cxx11=on.output new file mode 100644 index 0000000000..476923c248 --- /dev/null +++ b/tests/lac/linear_operator_13.with_cxx11=on.output @@ -0,0 +1,33 @@ + +DEAL::Block vector: v1: +DEAL::[block 0 ] 0.0 +DEAL::[block 1 ] 3.0 +DEAL::[block 2 ] 14. +DEAL::Block vector: v2: +DEAL::[block 0 ] 0.0 +DEAL::[block 1 ] 3.0 +DEAL::[block 2 ] 14. +DEAL::Block vector: v1: +DEAL::[block 0 ] 18. +DEAL::[block 1 ] 15. +DEAL::[block 2 ] 0.0 +DEAL::Block vector: v2: +DEAL::[block 0 ] 18. +DEAL::[block 1 ] 15. +DEAL::[block 2 ] 0.0 +DEAL::Block vector: v1: +DEAL::[block 0 ] 1.0 +DEAL::[block 1 ] 4.0 +DEAL::[block 2 ] 15. +DEAL::Block vector: v2: +DEAL::[block 0 ] 1.0 +DEAL::[block 1 ] 4.0 +DEAL::[block 2 ] 15. +DEAL::Block vector: v1: +DEAL::[block 0 ] 19. +DEAL::[block 1 ] 16. +DEAL::[block 2 ] 1.0 +DEAL::Block vector: v2: +DEAL::[block 0 ] 19. +DEAL::[block 1 ] 16. +DEAL::[block 2 ] 1.0 diff --git a/tests/lac/linear_operator_15.cc b/tests/lac/linear_operator_15.cc index 233b7cbdd7..18491944f4 100644 --- a/tests/lac/linear_operator_15.cc +++ b/tests/lac/linear_operator_15.cc @@ -49,9 +49,9 @@ int main() sparsity_pattern.copy_from(dsp); sparsity_pattern.compress(); - // | 2 1 0 | - // | 3 2 1 | - // | 4 3 2 | + // | 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) @@ -59,7 +59,7 @@ int main() auto op_a = linear_operator>(a); - auto triangulat_block_op = lower_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); + auto triangular_block_op = lower_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); BlockVector u; op_a.reinit_domain_vector(u, false); @@ -74,7 +74,7 @@ int main() v.block(i)[0] = 1; } - triangulat_block_op.vmult(v, u); + triangular_block_op.vmult(v, u); deallog << " -- vmult --" << std::endl; PRINTME("v", v); @@ -82,7 +82,7 @@ int main() for(unsigned int i = 0; i<3; ++i) v.block(i)[0] = 1; - triangulat_block_op.vmult_add(v, u); + triangular_block_op.vmult_add(v, u); deallog << " -- vmult_add --" << std::endl; PRINTME("v", v); @@ -90,7 +90,7 @@ int main() for(unsigned int i = 0; i<3; ++i) v.block(i)[0] = i+1; - triangulat_block_op.Tvmult(u, v); + triangular_block_op.Tvmult(u, v); deallog << " -- Tvmult --" << std::endl; PRINTME("u", u); @@ -98,7 +98,7 @@ int main() for(unsigned int i = 0; i<3; ++i) u.block(i)[0] = 1; - triangulat_block_op.Tvmult_add(u, v); + triangular_block_op.Tvmult_add(u, v); deallog << " -- Tvmult_add --" << std::endl; PRINTME("u", u); }