From 9a3a30c450754ce0bf276bec72ab714955c4528c Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Thu, 4 Jun 2015 23:40:20 +0200 Subject: [PATCH] rewrite old upper and triangular operators in terms of new ones --- include/deal.II/lac/block_linear_operator.h | 333 ++++++-------------- tests/lac/block_linear_operator_03.cc | 2 +- tests/lac/block_linear_operator_04.cc | 2 +- tests/lac/block_linear_operator_05.cc | 4 +- tests/lac/block_linear_operator_06.cc | 8 +- tests/lac/block_linear_operator_07.cc | 5 +- 6 files changed, 99 insertions(+), 255 deletions(-) diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h index 9118f12ce4..4d0e274fb4 100644 --- a/include/deal.II/lac/block_linear_operator.h +++ b/include/deal.II/lac/block_linear_operator.h @@ -325,8 +325,9 @@ block_diagonal_operator(const LinearOperator, - 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()) ); - - LinearOperator return_op; - - return_op.reinit_range_vector = [&block_matrix](Range &v, bool fast) - { - // Reinitialize the block vector to have the number of blocks - // equal to the number of row blocks of the matrix block_matrix. - v.reinit(block_matrix.n_block_rows()); - - // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < block_matrix.n_block_rows(); ++i) - linear_operator(block_matrix.block(i, 0)) - .reinit_range_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.reinit_domain_vector = [&block_matrix](Domain &v, bool fast) - { - // Reinitialize the block vector to have the number of blocks - // equal to the number of coloumn blocks of the matrix block_matrix. - v.reinit(block_matrix.n_block_cols()); - - // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) - linear_operator(block_matrix.block(0, i)) - .reinit_domain_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.vmult = [&block_matrix](Range &v, const Domain &u) - { - Assert(v.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); - Assert(u.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - - 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)); - for (unsigned int j = 1; j < i; ++j) - block_matrix.block(i,j).vmult_add(v.block(i), u.block(j)); - } - }; - - return_op.vmult_add = [&block_matrix](Range &v, const Domain &u) - { - Assert(v.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); - Assert(u.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - - for (unsigned int i = 1; i < block_matrix.n_block_rows(); ++i) - { - block_matrix.block(i,0).vmult_add(v.block(i), u.block(0)); - for (unsigned int j = 1; j < i; ++j) - block_matrix.block(i,j).vmult_add(v.block(i), u.block(j)); - } - }; - - return_op.Tvmult = [&block_matrix](Domain &v, const Range &u) - { - Assert(v.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); - Assert(u.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); - - for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) - { - 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)); - } - }; - - return_op.Tvmult_add = [&block_matrix](Domain &v, const Range &u) - { - Assert(v.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); - Assert(u.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); - - for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) - 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)); - }; - - return return_op; -} - - -/** - * @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 , @@ -481,7 +366,7 @@ lower_triangular_operator(const std::array | | a12 - * --------------- --------------- - * a20 | a21 | a22 | | - * @endcode + * This function is a specification of the above function that + * allows to work with block matrices @p block_matrix . * * @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 , +template , typename Domain = Range, typename BlockMatrix> LinearOperator -upper_triangular_operator(const BlockMatrix &); +lower_triangular_operator(const BlockMatrix &); -template +template LinearOperator -upper_triangular_operator(const BlockMatrix &block_matrix) +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()) ); - - LinearOperator return_op; - - return_op.reinit_range_vector = [&block_matrix](Range &v, bool fast) - { - // Reinitialize the block vector to have the number of blocks - // equal to the number of row blocks of the matrix block_matrix. - v.reinit(block_matrix.n_block_rows()); - - // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < block_matrix.n_block_rows(); ++i) - linear_operator(block_matrix.block(i, 0)) - .reinit_range_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.reinit_domain_vector = [&block_matrix](Domain &v, bool fast) - { - // Reinitialize the block vector to have the number of blocks - // equal to the number of coloumn blocks of the matrix block_matrix. - v.reinit(block_matrix.n_block_cols()); - - // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) - linear_operator(block_matrix.block(0, i)) - .reinit_domain_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.vmult = [&block_matrix](Range &v, const Domain &u) - { - Assert( v.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); - Assert( u.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - - for (unsigned int i = 0; i < block_matrix.n_block_rows() - 1; ++i) - { - block_matrix.block(i,block_matrix.n_block_rows() - 1).vmult(v.block(i), u.block(block_matrix.n_block_rows() - 1)); - for (unsigned int j = block_matrix.n_block_rows() - 2; j > i; --j) - block_matrix.block(i,j).vmult_add(v.block(i), u.block(j)); - } - v.block(block_matrix.n_block_rows() - 1) = 0; - }; - - return_op.vmult_add = [&block_matrix](Range &v, const Domain &u) - { - Assert( v.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); - Assert( u.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - - for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) - for (unsigned int j = i + 1; j < block_matrix.n_block_rows(); ++j) - block_matrix.block(i,j).Tvmult_add(v.block(i), u.block(j)); - }; - - return_op.Tvmult = [&block_matrix](Domain &v, const Range &u) - { - Assert( v.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); - Assert( u.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); - - for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) - { - 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)); - } - }; - - return_op.Tvmult_add = [&block_matrix](Domain &v, const Range &u) - { - Assert( v.n_blocks() == block_matrix.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); - Assert( u.n_blocks() == block_matrix.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); - - for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) - for (unsigned int j = 0; j < i; ++j) - block_matrix.block(j,i).vmult_add(v.block(i), u.block(j)); - }; + Assert(block_matrix.n_block_rows() == block_matrix.n_block_cols(), + ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) ); - return return_op; + 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 * - * This function is a specification of the above function that - * allows to work with std::array of std::array of 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> @@ -694,7 +500,7 @@ upper_triangular_operator(const std::array, + 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 * @@ -776,7 +618,8 @@ upper_triangular_operator(const std::array, +template , typename Domain = Range, typename BlockMatrix> LinearOperator @@ -785,33 +628,35 @@ block_triangular_inverse(const BlockMatrix &, bool lower = true); -template +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; + LinearOperator op_a; if (lower) { - op_a = lower_triangular_operator(block_matrix); + op_a = lower_triangular_operator(block_matrix); } else { - op_a = upper_triangular_operator(block_matrix); + 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)^block_matrix.n_block_cols() - for (unsigned int i = 0; i < block_matrix.n_block_cols() - 1; ++i) + // 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; diff --git a/tests/lac/block_linear_operator_03.cc b/tests/lac/block_linear_operator_03.cc index 21f89c8222..dad87eb770 100644 --- a/tests/lac/block_linear_operator_03.cc +++ b/tests/lac/block_linear_operator_03.cc @@ -59,7 +59,7 @@ int main() a.block(i,j).set(0, 0, 2 + i - j); auto op_a = linear_operator>(a); - auto triangular_block_op = lower_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); + auto triangular_block_op = lower_triangular_operator<3, BlockVector, BlockVector, BlockSparseMatrix >(a); BlockVector u; op_a.reinit_domain_vector(u, false); diff --git a/tests/lac/block_linear_operator_04.cc b/tests/lac/block_linear_operator_04.cc index e5bfcd2ed2..7b205309f5 100644 --- a/tests/lac/block_linear_operator_04.cc +++ b/tests/lac/block_linear_operator_04.cc @@ -59,7 +59,7 @@ int main() auto op_a = linear_operator>(a); - auto triangular_block_op = upper_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); + auto triangular_block_op = upper_triangular_operator<3, BlockVector, BlockVector, BlockSparseMatrix >(a); BlockVector u; op_a.reinit_domain_vector(u, false); diff --git a/tests/lac/block_linear_operator_05.cc b/tests/lac/block_linear_operator_05.cc index 09ed95701f..b4373fb389 100644 --- a/tests/lac/block_linear_operator_05.cc +++ b/tests/lac/block_linear_operator_05.cc @@ -57,8 +57,8 @@ int main() 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); + auto lower_triangular_block_op = lower_triangular_operator<3, BlockVector, BlockVector, BlockSparseMatrix >(a); + auto upper_triangular_block_op = upper_triangular_operator<3, BlockVector, BlockVector, BlockSparseMatrix >(a); BlockVector u1; diff --git a/tests/lac/block_linear_operator_06.cc b/tests/lac/block_linear_operator_06.cc index 9f46e5a045..69e9b65d5d 100644 --- a/tests/lac/block_linear_operator_06.cc +++ b/tests/lac/block_linear_operator_06.cc @@ -62,9 +62,9 @@ int main() d.set(i,i, 1.0 / (i+i +1) ); auto op_a = linear_operator< BlockVector >(a); - auto diagonal_inv = linear_operator< BlockVector >(d); + auto diagonal_inv = linear_operator< BlockVector, BlockVector >(d); // auto diagonal_inv = block_diagonal_operator(diag); - auto inverse_op_a = block_triangular_inverse< BlockVector, BlockVector, BlockSparseMatrix >(a, diagonal_inv); + auto inverse_op_a = block_triangular_inverse<3, BlockVector, BlockVector, BlockSparseMatrix >(a, diagonal_inv); auto identity = inverse_op_a * op_a; @@ -81,7 +81,7 @@ int main() u.block(i)[0] = 0;; v.block(i)[0] = 0; } - u.block(j)[0] = 1;; + u.block(j)[0] = 1; op_a.vmult(v, u); @@ -90,7 +90,7 @@ int main() deallog << " -- Inverse -- " << std::endl; inverse_op_a.reinit_domain_vector(u, false); - inverse_op_a.reinit_range_vector(v, 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) diff --git a/tests/lac/block_linear_operator_07.cc b/tests/lac/block_linear_operator_07.cc index d0c778628f..91c7ead21f 100644 --- a/tests/lac/block_linear_operator_07.cc +++ b/tests/lac/block_linear_operator_07.cc @@ -75,8 +75,8 @@ int main() 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); + auto lower_triangular_block_op = lower_triangular_operator<2, BlockVector, BlockVector, BlockSparseMatrix>(a); + auto upper_triangular_block_op = upper_triangular_operator<2, BlockVector, BlockVector, BlockSparseMatrix >(a); BlockVector u1; @@ -96,7 +96,6 @@ int main() PRINTME("v1", v1); PRINTME("v2", v2); - lower_triangular_block_op.vmult(v1, u1); upper_triangular_block_op.Tvmult(v2, u2); -- 2.39.5