From 399b54553b23dddcfe77a3d0fb0063d1f4e9691e Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Wed, 20 May 2015 18:52:01 +0200 Subject: [PATCH] preconditioner_from_diagonal_inverse preconditioner inverse check added --- include/deal.II/lac/linear_operator.h | 32 ++++ tests/lac/linear_operator_14.cc | 146 ++++++++++++++++++ .../linear_operator_14.with_cxx11=on.output | 0 3 files changed, 178 insertions(+) create mode 100644 tests/lac/linear_operator_14.cc create mode 100644 tests/lac/linear_operator_14.with_cxx11=on.output diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index c46a0a44e5..6cadd62cd2 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -2329,6 +2329,38 @@ upper_triangular_operator(BlockMatrix &a) return return_op; } + +template , + typename Domain = Range, + typename BlockMatrix> +LinearOperator +inverse_operator( BlockMatrix &a, + const LinearOperator &inverse_diagonal, + bool lower = true) +{ + + LinearOperator op_a; + + if(lower) + { + op_a = lower_triangular_operator(a); + } + else + { + op_a = upper_triangular_operator(a); + } + + auto id = identity_operator(op_a.reinit_range_vector); + auto result = identity_operator(op_a.reinit_range_vector); + + + for(unsigned int i = 0; i < a.n_block_cols() - 1; ++i) + result = id - inverse_diagonal * op_a * result; + + return result * inverse_diagonal; +} + + //@} DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/linear_operator_14.cc b/tests/lac/linear_operator_14.cc new file mode 100644 index 0000000000..a722ad1d7a --- /dev/null +++ b/tests/lac/linear_operator_14.cc @@ -0,0 +1,146 @@ +// --------------------------------------------------------------------- +// +// 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(12); + + // 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(); + + BlockSparseMatrix a (sparsity_pattern); + for (unsigned int i = 0; i < 3; ++i) + { + a.block(i,i).set(0, 0, i+i +1); + for (unsigned int j = 0; j < i; ++j) + a.block(i,j).set(0, 0, 10); + } + + BlockSparseMatrix d(sparsity_pattern); + for (unsigned int i = 0; i < 3; ++i) + d.set(i,i, 1.0 / (i+i +1) ); + // std::array diag {{a.block(0,0), a.block(1,1),a.block(2,2)}}; + + + auto op_a = linear_operator< BlockVector >(a); + auto diagonal_inv = linear_operator< BlockVector >(d); + // auto diagonal_inv = block_diagonal_operator(diag); + auto inverse_op_a = inverse_operator< BlockVector, BlockVector, BlockSparseMatrix >(a, diagonal_inv); + + auto identity = inverse_op_a * op_a; + + BlockVector u; + BlockVector v; + + deallog << " -- Matrix -- " << std::endl; + op_a.reinit_domain_vector(u, false); + op_a.reinit_range_vector(v, false); + for(unsigned int j = 0; j<3; ++j) + { + // check 1 + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0;; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1;; + + + op_a.vmult(v, u); + // identity.Tvmult(v2, u2); + + PRINTME("v", v); + } + + deallog << " -- Inverse -- " << std::endl; + inverse_op_a.reinit_domain_vector(u, false); + inverse_op_a.reinit_range_vector(v, false); + for(unsigned int j = 0; j<3; ++j) + { + // check 1 + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0;; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1;; + + + inverse_op_a.vmult(v, u); + // identity.Tvmult(v2, u2); + + PRINTME("v", v); + } + + deallog << " -- Identity -- " << std::endl; + identity.reinit_domain_vector(u, false); + identity.reinit_range_vector(v, false); + for(unsigned int j = 0; j<3; ++j) + { + // check 1 + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0;; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1;; + + + identity.vmult(v, u); + // identity.Tvmult(v2, u2); + + PRINTME("v", v); + } + + // op_a.vmult(v, u); + // // identity.Tvmult(v2, u2); + // + // PRINTME("v", v); + // + // inverse_op_a.vmult(v, u); + // // identity.Tvmult(v2, u2); + // + // PRINTME("v", v); + } +} diff --git a/tests/lac/linear_operator_14.with_cxx11=on.output b/tests/lac/linear_operator_14.with_cxx11=on.output new file mode 100644 index 0000000000..e69de29bb2 -- 2.39.5