From 22a90e08e511bcc60a12e367e790543b08918dd8 Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Thu, 4 Jun 2015 11:45:28 +0200 Subject: [PATCH] upper and lower triangular operator for array of array of linear operator --- include/deal.II/lac/block_linear_operator.h | 171 +++++++++++++++++++- tests/lac/block_linear_operator_07.cc | 10 -- 2 files changed, 168 insertions(+), 13 deletions(-) diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h index d17e290f33..9118f12ce4 100644 --- a/include/deal.II/lac/block_linear_operator.h +++ b/include/deal.II/lac/block_linear_operator.h @@ -399,7 +399,7 @@ lower_triangular_operator(const BlockMatrix &block_matrix) Assert(u.n_blocks() == block_matrix.n_block_cols(), ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - v.block(0) *= 0; + v.block(0) = 0; for (unsigned int i = 1; i < block_matrix.n_block_rows(); ++i) { block_matrix.block(i,0).vmult(v.block(i), u.block(0)); @@ -432,7 +432,7 @@ lower_triangular_operator(const BlockMatrix &block_matrix) for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) { - v.block(i) *= 0; + v.block(i) = 0; for (unsigned int j = i + 1; j < block_matrix.n_block_rows(); ++j) block_matrix.block(j,i).Tvmult_add(v.block(i), u.block(j)); } @@ -454,6 +454,90 @@ lower_triangular_operator(const BlockMatrix &block_matrix) } +/** + * @relates LinearOperator + * + * This function is a specification of the above function that + * allows to work with std::array of std::array of LinearOperator + * + * @ingroup LAOperators + */ + +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][0].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[0][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 * @@ -561,7 +645,7 @@ upper_triangular_operator(const BlockMatrix &block_matrix) for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) { - v.block(i) *= 0; + v.block(i) = 0; for (unsigned int j = 0; j < i; ++j) block_matrix.block(j,i).vmult_add(v.block(i), u.block(j)); } @@ -583,6 +667,87 @@ upper_triangular_operator(const BlockMatrix &block_matrix) } +/** + * @relates LinearOperator + * + * This function is a specification of the above function that + * allows to work with std::array of std::array of LinearOperator + * + * @ingroup LAOperators + */ + +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][0].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[0][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 * diff --git a/tests/lac/block_linear_operator_07.cc b/tests/lac/block_linear_operator_07.cc index d902ad86d5..d0c778628f 100644 --- a/tests/lac/block_linear_operator_07.cc +++ b/tests/lac/block_linear_operator_07.cc @@ -90,16 +90,6 @@ int main() // check 1 - // for(unsigned int i = 0; i<2; ++i) - // { - // for(unsigned int j = 0; j<2; ++j) - // { - // u1.block(i)[j] = i+1; - // u2.block(i)[j] = i+1; - // v1.block(i)[j] = 1; - // v2.block(i)[j] = 1; - // } - // } reinit_vec(u1,u2,v1,v2); PRINTME("u1", u1); PRINTME("u2", u2); -- 2.39.5