From: ESeNonFossiIo Date: Tue, 16 Jun 2015 10:56:45 +0000 (+0200) Subject: reneaming and improving documentation X-Git-Tag: v8.3.0-rc1~92^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e03a4ae264e6ac33e59c8c2c2990fa27505d5436;p=dealii.git reneaming and improving documentation --- diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h index 3ed5833ab2..1a480a7d45 100644 --- a/include/deal.II/lac/block_linear_operator.h +++ b/include/deal.II/lac/block_linear_operator.h @@ -665,21 +665,36 @@ block_triangular_inverse(const BlockMatrix &block_matrix, /** * @relates LinearOperator * - * This function implement a gauss elimination argument to invert a lower + * 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 * @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 + * matrix: + * A00 x0 + ... = y0 + * A01 x0 + A11 x1 + ... = y1 + * ... ... + * A0n x0 + A1n x1 + ... + Ann xn = yn + * + * First of all, x0 = A00^-1 y0 + * Then, we can use x0 to recover x1: + * x1 = A11^-1 ( y1 - A01 x0 ) + * and therefore: + * xn = Ann^-1 ( yn - A0n x0 - ... - A(n-1)n x(n-1) ) + * * Caveat: Tvmult and Tvmult_add have not been implemented, yet. This may lead to mistakes. * @ingroup LAOperators */ -template -LinearOperator -gauss_lower_block_inverse(const BlockMatrix &block_matrix, - const std::array, n> &inverse_diagonal) +template +LinearOperator +block_forward_substitution(const BlockMatrix &block_matrix, + const std::array, n> &inverse_diagonal) { - LinearOperator return_op; + LinearOperator return_op; return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast) { @@ -693,7 +708,7 @@ gauss_lower_block_inverse(const BlockMatrix &block_matrix, v.collect_sizes(); }; - return_op.reinit_domain_vector = [inverse_diagonal](Domain &v, bool fast) + return_op.reinit_domain_vector = [inverse_diagonal](Range &v, bool fast) { // Reinitialize the block vector to m blocks: v.reinit(n); @@ -705,7 +720,7 @@ gauss_lower_block_inverse(const BlockMatrix &block_matrix, v.collect_sizes(); }; - return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Domain &u) + return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Range &u) { static GrowingVectorMemory vector_memory; typename Range::BlockType *tmp = vector_memory.alloc(); @@ -725,7 +740,7 @@ gauss_lower_block_inverse(const BlockMatrix &block_matrix, vector_memory.free(tmp); }; - return_op.vmult_add = [&block_matrix, inverse_diagonal](Range &v, const Domain &u) + return_op.vmult_add = [&block_matrix, inverse_diagonal](Range &v, const Range &u) { static GrowingVectorMemory vector_memory; typename Range::BlockType *tmp = vector_memory.alloc(); @@ -751,21 +766,37 @@ gauss_lower_block_inverse(const BlockMatrix &block_matrix, /** * @relates LinearOperator * - * This function implement a gauss elimination argument to invert an upper + * This function implement a back substitution argument to invert an upper * block triangular matrix. * It takes as argement a block matrix @p 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 + * matrix: + + * A00 x0 + A01 x1 + ... + A0n xn = yn + * A11 x1 + ... = y1 + * ... .. + * Ann xn = yn + * + * First of all, xn = Ann^-1 yn + * Then, we can use xn to recover x(n-1): + * x(n-1) = A(n-1)(n-1)^-1 ( y(n-1) - A(n-1)n x(n-1) ) + * and therefore: + * x0 = A00^-1 ( y0 - A0n xn - ... - A01 x1 ) + * * Caveat: Tvmult and Tvmult_add have not been implemented, yet. This may lead to mistakes. * @ingroup LAOperators */ -template -LinearOperator -gauss_upper_block_inverse(const BlockMatrix &block_matrix, - const std::array, n> &inverse_diagonal) +template +LinearOperator +block_back_substitution(const BlockMatrix &block_matrix, + const std::array, n> &inverse_diagonal) { - LinearOperator return_op; + LinearOperator return_op; return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast) { @@ -779,7 +810,7 @@ gauss_upper_block_inverse(const BlockMatrix &block_matrix, v.collect_sizes(); }; - return_op.reinit_domain_vector = [inverse_diagonal](Domain &v, bool fast) + return_op.reinit_domain_vector = [inverse_diagonal](Range &v, bool fast) { // Reinitialize the block vector to m blocks: v.reinit(n); @@ -791,7 +822,7 @@ gauss_upper_block_inverse(const BlockMatrix &block_matrix, v.collect_sizes(); }; - return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Domain &u) + return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Range &u) { static GrowingVectorMemory vector_memory; typename Range::BlockType *tmp = vector_memory.alloc(); @@ -811,7 +842,7 @@ gauss_upper_block_inverse(const BlockMatrix &block_matrix, vector_memory.free(tmp); }; - return_op.vmult_add = [&block_matrix, inverse_diagonal](Range &v, const Domain &u) + return_op.vmult_add = [&block_matrix, inverse_diagonal](Range &v, const Range &u) { static GrowingVectorMemory vector_memory; typename Range::BlockType *tmp = vector_memory.alloc(); diff --git a/tests/lac/block_linear_operator_09.cc b/tests/lac/block_linear_operator_09.cc index 49c726f56f..ad4532af79 100644 --- a/tests/lac/block_linear_operator_09.cc +++ b/tests/lac/block_linear_operator_09.cc @@ -13,7 +13,7 @@ // // --------------------------------------------------------------------- -// Test gauss_upper_block_inverse and gauss_lower_block_inverse: +// Test block_back_substitution and block_forward_substitution: #include "../tests.h" @@ -65,9 +65,8 @@ int main() 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 diagonal_inv = block_diagonal_operator(diag); - auto inverse_op_a = gauss_lower_block_inverse<3, BlockVector, BlockVector, BlockSparseMatrix >(a, {{d00, d11, d22}}); - // + + auto inverse_op_a = block_forward_substitution< BlockSparseMatrix >(a, {{d00, d11, d22}}); auto identity = inverse_op_a * op_a; BlockVector u; @@ -153,9 +152,8 @@ int main() 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 diagonal_inv = block_diagonal_operator(diag); - auto inverse_op_a = gauss_upper_block_inverse<3, BlockVector, BlockVector, BlockSparseMatrix >(a, {{d00, d11, d22}}); - // + + auto inverse_op_a = block_back_substitution< BlockSparseMatrix >(a, {{d00, d11, d22}}); auto identity = inverse_op_a * op_a; BlockVector u;