From: ESeNonFossiIo Date: Wed, 20 May 2015 08:18:44 +0000 (+0200) Subject: lower triangular block matrix X-Git-Tag: v8.3.0-rc1~155^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9e4c4d88872b08b3bf6bc59bcc52c4627167c8cc;p=dealii.git lower triangular block matrix lower triangular matrix tests --- diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index d02f153fe8..0f39f75138 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -2126,6 +2126,109 @@ operator*(const PackagedOperation &comp, return return_comp; } + + +template , + typename Domain = Range, + typename BlockMatrix> +LinearOperator +lower_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())); + + v.block(0) *= 0; + for (unsigned int i = 1; i < a.n_block_rows(); ++i) + { + a.block(i,0).vmult(v.block(i), u.block(0)); + for (unsigned int j = 1; j < i; ++j) + a.block(i,j).vmult_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 = 1; i < a.n_block_rows(); ++i) + { + a.block(i,0).vmult_add(v.block(i), u.block(0)); + for (unsigned int j = 1; j < i; ++j) + a.block(i,j).vmult_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 = i + 1; j < a.n_block_rows(); ++j) + a.block(j,i).Tvmult_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 = i + 1; j < a.n_block_rows(); ++j) + a.block(j,i).Tvmult_add(v.block(i), u.block(j)); + }; + + return return_op; +} + + //@} DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/linear_operator_11.with_cxx11=on.output b/tests/lac/linear_operator_11.with_cxx11=on.output new file mode 100644 index 0000000000..1507567775 --- /dev/null +++ b/tests/lac/linear_operator_11.with_cxx11=on.output @@ -0,0 +1,21 @@ + +DEAL:: -- vmult -- +DEAL::Block vector: v: +DEAL::[block 0 ] 0.0 +DEAL::[block 1 ] 3.0 +DEAL::[block 2 ] 10. +DEAL:: -- vmult_add -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.0 +DEAL::[block 1 ] 4.0 +DEAL::[block 2 ] 11. +DEAL:: -- Tvmult -- +DEAL::Block vector: u: +DEAL::[block 0 ] 18. +DEAL::[block 1 ] 9.0 +DEAL::[block 2 ] 0.0 +DEAL:: -- Tvmult_add -- +DEAL::Block vector: u: +DEAL::[block 0 ] 19. +DEAL::[block 1 ] 10. +DEAL::[block 2 ] 1.0 diff --git a/tests/lac/linear_operator_15.cc b/tests/lac/linear_operator_15.cc new file mode 100644 index 0000000000..233b7cbdd7 --- /dev/null +++ b/tests/lac/linear_operator_15.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 1 0 | + // | 3 2 1 | + // | 4 3 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 + i - j); + + auto op_a = linear_operator>(a); + + auto triangulat_block_op = lower_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; + } + + triangulat_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; + + triangulat_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; + + triangulat_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; + + triangulat_block_op.Tvmult_add(u, v); + deallog << " -- Tvmult_add --" << std::endl; + PRINTME("u", u); + } +}