From: ESeNonFossiIo Date: Mon, 15 Jun 2015 22:04:39 +0000 (+0200) Subject: gauss elimination argument for linear operator X-Git-Tag: v8.3.0-rc1~92^2~7 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c41586bed57f33e82e110e3e10683a3f159032a0;p=dealii.git gauss elimination argument for linear operator --- diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h index 4d0e274fb4..3ed5833ab2 100644 --- a/include/deal.II/lac/block_linear_operator.h +++ b/include/deal.II/lac/block_linear_operator.h @@ -662,6 +662,177 @@ block_triangular_inverse(const BlockMatrix &block_matrix, return result * inverse_diagonal; } +/** + * @relates LinearOperator + * + * This function implement a gauss elimination 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. + * + * 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) +{ + LinearOperator return_op; + + return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(n); + + // And reinitialize every individual block with reinit_range_vectors: + for (unsigned int i = 0; i < n; ++i) + inverse_diagonal[i].reinit_range_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.reinit_domain_vector = [inverse_diagonal](Domain &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(n); + + // And reinitialize every individual block with reinit_domain_vectors: + for (unsigned int i = 0; i < n; ++i) + inverse_diagonal[i].reinit_domain_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Domain &u) + { + static GrowingVectorMemory vector_memory; + typename Range::BlockType *tmp = vector_memory.alloc(); + + inverse_diagonal[0].vmult(v.block(0), u.block(0)); + for (unsigned int i=1; i vector_memory; + typename Range::BlockType *tmp = vector_memory.alloc(); + + inverse_diagonal[0].vmult_add(v.block(0), u.block(0)); + for (unsigned int i=1; i +LinearOperator +gauss_upper_block_inverse(const BlockMatrix &block_matrix, + const std::array, n> &inverse_diagonal) +{ + LinearOperator return_op; + + return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(n); + + // And reinitialize every individual block with reinit_range_vectors: + for (unsigned int i = 0; i < n; ++i) + inverse_diagonal[i].reinit_range_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.reinit_domain_vector = [inverse_diagonal](Domain &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(n); + + // And reinitialize every individual block with reinit_domain_vectors: + for (unsigned int i = 0; i < n; ++i) + inverse_diagonal[i].reinit_domain_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.vmult = [&block_matrix, inverse_diagonal](Range &v, const Domain &u) + { + static GrowingVectorMemory vector_memory; + typename Range::BlockType *tmp = vector_memory.alloc(); + + inverse_diagonal[n-1].vmult(v.block(n-1),u.block(n-1)); + for (int i=n-2; i>=0; --i) + { + inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true); + *tmp = u.block(i); + *tmp *= -1.; + for (int j=i+1; j vector_memory; + typename Range::BlockType *tmp = vector_memory.alloc(); + + inverse_diagonal[n-1].vmult_add(v.block(n-1),u.block(n-1)); + for (int i=n-2; i>=0; --i) + { + inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true); + *tmp = u.block(i); + *tmp *= -1.; + for (int j=i+1; j